y-  f5'-c-6q  *  t 

The  Pennsylvania  State  University 
APPLIED  RESEARCH  LABORATORY 
P.  0.  Box  30 

State  College,  PA  16804 


A  COUPLED,  PARABOLIC -MARCHING  METHOD  FOR  THE 
PREDICTION  OF  THREE-DIMENSIONAL  VISCOUS 
INCOMPRESSIBLE  TURBOHACHINERY  FLOWS 


by 

Kevin  Richard  Kircley 


$ 


DT! 

elec 

0CT  2  7 1988 


y 


<V 


Technical  Report  No.  TR  88-014 
October  1988 


Supported  by: 

Naval  Sea  Systems  Command 


L.  R.  Hettche,  Director 
Applied  Research  Laboratory 


Approved  for  public  release;  distribution  unlimited 


i 


w:iia<nw»> 


U  REPORT  SECURITY  CLASSIFICATION 

Unclassified 


2a.  SECURITY  CLASSIFICATION  authority 


2b  OECLASSIFICATION/OOWNGRAOING  schedule 


4  PERFORMING  ORGANIZATION  REPORT  NUM8£R(Si 


REPORT  DOCUMENTATION  PAGE 


lb.  RESTRICTIVE  MARKINGS 


3.  DISTRIBUTION/ AVAILA8IUTY  OF  REPORT 


Unlimited 


S.  MONITORING  ORGANIZATION  REPORT  NUM8ER(S) 


x  ADORESS  (Cry.  Slat*.  and  ZIP  Code) 

P.  0.  Box  30 

State  College,  PA  16804 


la.  NAME  of  funding  /sponsoring 

ORGANIZATION 

Naval  Sea  Systems  Command 


ic.  ADDRESS  fOTy.  Stitt,  and  ZIP  Cod*) 
Department  of  the  Navy 
Washington,  DC  20362 


6b  OFFICE  SYMBOL  7*.  NAME  OF  MONITORING  ORGANIZATION 

(if  applicant)  David  Taylor  Research  Center 

y  ARL  Department  of  the  Naw 


7b.  ADDRESS  (City.  Stitt,  and  ZIP  Codt) 

Eethesda,  HD  20084 


Sb.  OFFICE  SYMBOL  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

■KiSSF*" 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM  PROJECT 

ELEMENT  NO.  NO.  • 


i  title  (Inciudt  Stcurrty  damnation) 

A  Coupled,  Parabolic-Marching  Method  for  the  Prediction  of  Three-Dimensional 
Viscous  Incompressible  Turbomachinerv  Flows  (Unclassified) 


-2  PERSONAL  AUTHOR(S) 

Kevin  Richard  Kirtlev 


Ua  type  OF  REPORT 

Ph.  D.  Thesis 


L  6  SUPPLEMENTARY  NOTATION 


COSATI  COOES 


13b.  TIME  COVERED 
FROM _  TO 


14  DATF  OF  REPORT  (Ytir.  Month.  Day) 

October  1988 


9  ABSTRACT  ( Continut  on  reverie  if  ntetoary  and  idtntity  by  block  number) 

— ^  A  new  coupled  parabolic-marching  method  was  developed  to  solve  the 
three-dimensional  incompressible  Navier-Stokes  equation  for  turbulent 
turbomachinery  flows.  Earlier  space-marching  methods  were  analyzed  to 
determine  their  global  stability  during  multiple  passes  of  the 
computational  domain.  The  methods  were  found  to  be  unconditionally 

t 

unstable  even  when  an  extra  equation  for  the  pressure,  namely  the 
Poisson  equation  for  the  pressure,  was  used  between  passes  of  the 


70  C"i  <  nlbuiiON  .'AVAILABILITY  OF  ABSTRACT  21.  ABSTRACT  SECURITY  CLASSIFICATION 

IS unciassifiedajnlimited  □  same  as  RPT.  □  otic  users  Unclassified 


a  NAME  OF  RESPONSIBLE  INOIVIOUAi  122b.  TELEPHONE  (JnOude  Area  Coda)  I  22c  OFFICE  SYMBOL 


O  FORM  1473, 84  MAR 


83  APR  tdition  may  b«  mad  until  exhausted. 
All  other  editions  are  obsolete. 


Unclassified 


ii 


Unclassified 


ueumrr  claaiificatio*  or  thii  r*«c 


domain.  Relaxation  of  one  constraint  curing  tr.e  solution  process  was 
founo  tc  be  necessary  for  the  successful  calculation  of  a  complex  flow. 
Thus,  the  method  of  pseudocompressibiiity  was  i.-.troouceo  into  the 
partially  parabolized  Navier-Stokes  equation  to  relax  the  mass  flow 
constraint  during  the  iterative  process.  This  new  method  was  found  to 
be  stable  during  a  forward-marching  integration  as  well  as  globally 
stable  during  successive  passes  of  the  domain.  With  consistent 

discretization,  the  new  method  was  found  to  be  convergent. 

Also  investigated  was  the  splitting  error  which  arises  from  the  use 

of  the  LBI  scheme  in  a  three-dimensional  parabolic-marching  method.  The 
splitting  error  was  found  to  be  extremely  important  for  coarse  grid 
computations  and  was  analyzed  and  demonstrated  for  a  strongly  curved 
duct  flow.  A  simple  iterative  solution  method  was  developed  which 
reduces  the  effect  of  the  splitting  error  for  three-dimensional 
computations. 

Several  turbulence  models  were  introduced  for  the  computation  of 
turbulent  flows.  The  three  models  used  were  the  algebraic  eddy 
viscosity  model,  the  two  equation  k-e  model,  and  the  algebraic  Reynolds 
stress  model  (ARSM)  which  introduces  the  strong  effect  of  rotation  on 
the  turbulence  structure. 

The  new  parabolic -marching  method  using  the  new  solution  method  was 
used  to  compute  several  flows  including  the  laminar  and  turbulent  flow 
in  an  S-shaped  duct  and  the  turbulent  flow  in  an  end  wall  cascade  and 
compressor  rotor.  For  the  laminar  flow,  the  agreement  with  the  analysis 
and  the  experimental  data  was  excellent.  For  the  turbulent  flows,  the 
pressu-e  distributions  were  accurately  predicted  as  were  blade  and  end 
wall  boundary  layers  and  wakes.  Prediction  accuracy  for  the  rotor  flow 
was  adequate  with  good  resolution  of  the  suction  side  boundary  layers, 
secondary  flow,  and  pressure  losses  at  the  rotor  exit. 


Unclassified 


SCCUMTV  CLAUinCATION  or  THIS  PAOK 


■ 


I 


I 


ABSTRACT 


iii 


A  new  coupled  parabolic-marching  method  was  developed  to  solve  the 
three-dimensional  incompressible  Navier-Stokes  equation  for  turbulent 
turbomachinery  flows.  Earlier  space-marching  methods  were  analyzed  to 
determine  their  global  stability  during  multiple  passes  of  the 
computational  domain.  The  methods  were  found  to  be  unconditionally 
unstable  even  when  an  extra  equation  for  the  pressure,  namely  the 
Poisson  equation  for  the  pressure,  was  used  between  passes  of  the 
domain.  Relaxation  of  one  constraint  during  the  solution  process  was 
found  to  be  necessary  for  the  successful  calculation  of  a  complex  flow. 
Thus,  the  method  of  pseudocompressibility  was  introduced  into  the 
partially  parabolized  Navier-Stokes  equation  to  relax  the  mass  flow 
constraint  during  the  iterative  process.  This  new  method  was  found  to 
be  stable  during  a  forward-marching  integration  as  well  as  globally 
stable  during  successive  passes  of  the  domain.  With  consistent 
discretization,  the  new  method  was  found  to  be  convergent. 

Also  investigated  was  the  splitting  error  which  arises  from  the  use 
of  the  LBI  scheme  in  a  three-dimensional  parabolic-marching  method.  The 
splitting  error  was  found  to  be  extremely  important  for  coarse  grid 
computations  and  was  analyzed  and  demonstrated  for  a  strongly  curved 
duct  flow.  A  simple  iterative  solution  method  was  developed  which 
reduces  the  effect  of  the  splitting  error  for  three-dimensional 
computations. 

Several  turbulence  models  were  introduced  for  the  computation  of 


i 


turbulent  flows.  The  three  models  used  were  the  algebraic  eddy 
viscosity  model,  the  two  equation  k-e  model,  and  the  algebraic  Reynolds 
stress  model  (ARSM)  which  introduces  the  strong  effect  of  rotation  on 


the  turbulence  structure. 

The  new  parabolic-marching  method  using  the  new  solution  method  was 
used  to  compute  several  flows  including  the  laminar  and  turbulent  flow 
in  an  S-shaped  duct  and  the  turbulent  flow  in  an  end  wall  cascade  and 
compressor  rotor.  For  the  laminar  flow,  the  agreement  with  the  analysis 
and  the  experimental  data  was  excellent.  For  the  turbulent  flows,  the 
pressure  distributions  were  accurately  predicted  as  were  blade  and  end 
wall  boundary  layers  and  wakes.  Prediction  accuracy  for  the  rotor  flow 
was  adequate  with  good  resolution  of  the  suction  side  boundary  layers, 
secondary  flow,  and  pressure  losses  at  the  rotor  exit. 


|  Accession  For  j 

HTIS  GRA4I 

DTIC  TAB 

fa 

Unannounced 

□ 

Justification _ 

V 


TABLE  OF  CONTENTS 

page 

ABSTRACT  .  . .  iii 

LIST  OF  FIGURES .  vii 

LIST  OF  SYMBOLS . xii 

ACKNOWLEDGMENTS  .  xv 

CHAPTER 

1  INTRODUCTION  .  1 

1 . 1  The  Problem .  1 

1.2  Review  of  Related  Studies  .  3 

1.2.1  Time-Marching  Methods  .  3 

1.2.2  Space-Marching  Methods  .  6 

1.2.3  Parabolic-Marching  Methods  .  10 

1.2.4  Elliptic  Methods  .  12 

1.3  Objectives  and  Method  of  Approach .  13 

2  ANALYSIS  AND  MODIFICATION  OF  EXISTING  TECHNIQUES  .  16 

2.1  Introduction .  16 

2.2  Difficulties  With  Existing  Techniques  ....  17 

2.3  Method  of  Pouagare  and  Lakshminarayana  (PL)  .  19 

2.3.1  Eigenvalue  Analysis  .  19 

2.3.2  Global  Stability  Analysis  .  22 

2.3.3  Remarks  on  the  Method .  31 

2.4  Modified  PL  Method  (MPL) .  33 

2.4.1  Eigenvalue  Analysis  .  34 

2.4.2  Global  Stability  Analysis  .  36 

2.4.3  Remarks  on  the  Method .  40 

2.5  Poisson  Equation  for  the  Pressure .  43 

2.6  Conclusion .  53 

3  NEW  PARABOLIC-MARCHING  METHOD  (NPM)  .  55 

3.1  Introduction .  55 

3.2  Eigenvalue  Analysis  .  57 

3.3  Global  Stability  Analysis  .  59 

3.4  Convergence  of  the  Method .  68 

3.5  Transformation  and  Discretization  of  the 

Equation .  72 

3.6  New  Solution  Procedure .  75 


vi 


3.6.1  Splitting  Error  in  LBI  Scheme .  75 

3.6.2  New  Solution  Method .  78 

3.6.3  Error  in  Curved  Duct  Flow .  80 

3.7  Boundary  Conditions  .  88 

3.8  Algorithm  Verification  .  90 

3.8.1  Developing  Laminar  Flow  in  a  Straight  Duct  91 

3.8.2  Laminar  Flow  in  an  S-Shaped  Duct  ....  97 

3.9  Closing  Remarks .  107 

4  TURBULENCE  CLOSURE  SCHEMES  .  Ill 

4.1  Introduction .  Ill 

4.2  Algebraic  Eddy  Viscosity  Model  .  112 

4.3  Two-Equation  Model  .  115 

4.3.1  Boundary  Conditions  .  118 

4.3.2  Initial  Conditions  .  119 

4.4  Algebraic  Reynolds  Stress  Model  .  121 

4.5  Methods  of  Solution .  123 

5  COMPUTATION  OF  COMPLEX  TURBOMACHINERY  FLOWS  ...  124 

5.1  Introduction .  124 

5.2  Flow  in  an  S-Shaped  Duct .  124 

5.3  Flow  in  an  End  Wall  Cascade .  136 

5.4  Compressor  Cascade  Wake  Flow .  152 

5.5  Flow  in  an  Axial  Flow  Compressor  Rotor  ....  163 

6  CONCLUDING  REMARKS  .  186 

6.1  Summary  of  Conclusions .  186 

6.2  Suggestions  for  Future  Research  .  1Q0 

REFERENCES .  191 


196 


APPENDIX:  COMPUTATIONAL  GRID  GENERATION 


vii 


LIST  OF  FIGURES 


page 


1 . 1  Nature  of  the  Flow  in  a  Turbomachine .  2 

2.1  Real  Part  of  X*X,  o=-0.01  for  All  Wave  Numbers 

for  the  PL  method  (u=v=1 , Ax=Ay=1 ,  Re=1000)  .  27 

2.2  Real  Part  of  X*X,  o=-1.00  for  All  Wave  Numbers 

for  the  PL  method  (u=v=1 , Ax=Ay=1 ,  Re=1000)  .  28 

2.3  Real  Part  of  X*X,  o=-0.01  for  All  Wave  Numbers 

for  the  PL  method  (u=v=1 , Ax=0.01  Ay=1,  Re=1000)  .  .  29 

2.4  Real  Part  of  X*X,  o=-0.01  for  All  Wave  Numbers 

for  the  PL  method  (u=v=1 ,Ax=100  Ay= 1 ,  Re=1000)  ...  30 

2.5  Real  Part  of  X*X,  b=-100  for  All  Wave  Numbers 

for  the  MPL  method  (u=v=1 ,Ax=Ay=1 ,  Re=1000)  ....  38 

2.6  Real  Part  of  X*X,  b=-100  for  All  Wave  Numbers 

for  the  MPL  method  (u=v=1 ,Ax=0.01  Ay=1,  Re=1000)  .  .  39 

2.7  Real  Part  of  X*X,  b=-100  for  All  Wave  Numbers 

for  the  MPL  method  (u=v=1 ,Ax=100  Ay=1,  Re=1000)  .  .  41 

2.8  Convergence  History  for  the  MPL  method  for  Developing 

Flow  in  a  Square  Duct,  Re=100 .  42 

2.9  Geometry  of  the  S-Shaped  Duct .  45 

2.10  Streamwise  Velocity  Profiles  for  the  S-Duct, 

Station  2,  MPL  method,  Re=40000  .  46 

2.11  Transverse  Velocity  Profiles  for  the  S-Duct, 

Station  2,  MPL  method,  Re=40000  .  47 

2.12  Streamwise  Velocity  Profiles  for  the  S-Duct, 

Station  5,  MPL  method,  Re=40000  .  48 

2.13  Transverse  Velocity  Profiles  for  the  S-Duct, 

Station  5,  MPL  method,  Re=40000  .  49 

2.14  Pressure  Computed  from  the  Poisson  Equation  Using 

the  MPL  method,  Re=40000  .  51 

3.1  Real  Part  of  X*Xi,  a=0.2,  for  All  Wave  Numbers 

for  the  NPM  method  (u=v=1 , Ax=Ay=1 ,  Re=1000)  ....  62 

3.2  Real  Part  of  X*Xi,  a=0.2,  for  All  Wave  Numbers 

for  the  NPM  method  (u=v=1 , Ax=0.01  Ay=1,  Re=1000)  .  .  64 


viii 


3.3  Real  Part  of  X*Xi,  a=0.2,  for  All  Wave  Numbers 

for  the  NPM  method  (u=v=1 , Ax=100  Ay=1,  Re=1000)  .  .  65 

3.4  Real  Part  of  X*X -j  for  Various  Values  of  a,  for  All 
Wave  Numbers  for  the  NPM  method  (u=v=1 , Ax=Ay=1 , 

Re=1000) . .  66 

3.5  Convergence  History  to  Machine  Precision  .  69 

3.6  Real  Part  of  X*Xi,  a=0.5,  for  All  Wave  Numbers 

for  the  NPM  method  (u=v=1  ,Ax=Ay=1 ,  Re=1000)  ....  70 

3.7  Strongly  Curved  Duct  Geometry  .  81 

3.8  Fine  Grid  Results  of  Pouagare  and  Lakshminarayana 

(1986)  for  the  Strongly  Curved  Duct  Flow .  82 

3.9  Computed  Streamwise  Velocity  Profiles  for  the 

Strongly  Curved  Duct  Flow,  PL  method  .  83 

3.10  Computed  Transverse  Velocity  Profiles  for  the 

Strongly  Curved  Duct  Flow,  PL  method .  84 

3.11  Computed  Residual  in  the  Continuity  Equation  for  the 

Three  Solution  Methods  .  86 

3.12  Convergence  History  for  the  Iterative  LBI  Scheme  .  87 

3.13  Convergence  History  for  the  Straight  Duct .  92 

3.14  Computed  Centerline  Coefficient  of  Pressure  for 

the  Straight  Duct .  94 

3.15  Computed  Centerline  Velocity  for  the  Straight  Duct  95 

3.16  Computed  Fully  Developed  Velocity  Profiles  for 

the  Straight  Duct .  96 

3.17  Convergence  History  for  the  S-Duct .  98 

3.18  Computed  Coefficient  of  Pressure  for  the  S-Duct  .  .  100 

3.19  Streamwise  Velocity  Profiles  for  the  S-Duct  at 

Station  1 .  101 

3.20  Transverse  Velocity  Profiles  for  the  S-Duct  at 

Station  1 .  102 

3.21  Streamwise  Velocity  Profiles  for  the  S-Duct  at 

Station  2 .  103 

3.22  Transverse  Velocity  Profiles  for  the  S-Duct  at 

Station  2 .  104 


ix 


3.23  Streamwise  Velocity  Profiles  for  the  S-Duct  at 

Station  4 .  105 

3.24  Transverse  Velocity  Profiles  for  the  S-Duct  at 

Station  4 .  106 

3.25  Streamwise  Velocity  Profiles  for  the  S-Duct  at 

Station  5 .  108 

3.26  Transverse  Velocity  Profiles  for  the  S-Duct  at 

Station  5 .  109 


5.1  Convergence  History  for  the  Turbulent  S-Duct  ....  126 

5.2  Streamwise  Velocity  Profiles  for  the  Turbulent  S-Duct 


at  Station  1 .  127 

5.3  Transverse  Velocity  Profiles  for  the  Turbulent  S-Duct 

at  Station  1 .  128 

5.4  Streamwise  Velocity  Profiles  for  the  Turbulent  S-Duct 

at  Station  2 .  129 

5.5  Transverse  Velocity  Profiles  for  the  Turbulent  S-Duct 

at  Station  2 .  13 1 

5.6  Streamwise  Velocity  Profiles  for  the  Turbulent  S-Duct 

at  Station  4 .  132 

5.7  Transverse  Velocity  Profiles  for  the  Turbulent  S-Duct 

at  Station  4 .  1 33 

5-8  Streamwise  Velocity  Profiles  for  the  Turbulent  S-Duct 

at  Station  5 .  1 34 

5.9  Transverse  Velocity  Profiles  for  the  Turbulent  S-Duct 

at  Station  5 .  1 35 

5.10  Computed  Coefficient  of  Pressure  for  the  Turbulent 

S-Duct .  137 


5.11  Geometry  and  Grid  for  the  End  Wall  Cascade  ....  1 38 

5.12  Convergence  History  for'  the  End  Wall  Cascade  .  .  .  140 

5.13  Pressure  Distribution  for  the  End  Wall  Cascade  .  .  1 4 1 


5.14  Streamwise  Velocity  Profiles  at  44%  Chord  for  the 

tnd  Wall  Cascade .  142 

5.15  Streamwise  Velocity  Profiles  at  88%  Chord  for  the 

End  Wall  Cascade .  143 


X 


5.16  Streamwise  Momentum  Thickness  for  the  End  Wall 

Cascade  . .  145 

5.17  Transverse  Velocity  Profiles  at  44?  Chord  for  the 

End  Wall  Cascade .  146 

5.18  Transverse  Velocity  Profiles  at  88?  Chord  for  the 

End  Wall  Cascade .  147 

5.19  Secondary  Flow  Velocities  at  44?  Chord  for  the  End 

Wall  Cascade .  148 

5.20  Secondary  Flow  Velocities  at  66?  Chord  for  the  End 

Wall  Cascade .  149 

5.21  Secondary  Flow  Velocities  at  88?  Chord  for  the  End 

Wall  Cacsade .  150 

5.22  Computed  Mass  A  eraged  Secondary  Flow  Kinetic  Energy 

for  the  End  Wall  Cascade .  151 

5.23  Geometry  and  Grid  for  Hobbs  Cascade  Wake  Flow  ...  153 

5.24  Convergence  History  f<v  Hobbs  Cascade,  Alg.  Model  .  154 

5-25  Centerline  Velocity  Distribution  for  Hobbs  Cascade, 

Alg.  Model .  155 


5.26  Wake  Profiles  for  Hobbs  Cascade,  Alg.  Model  ....  157 

5.27  Convergence  History  for  Hobbs  Cascade,  k-e  Model  .  158 

5.28  Centerline  Velocity  Distribution  for  Hobbs  Cacsade, 


k-e  Model .  159 

5.29  Wake  Profiles  for  Hobbs  Cacsade,  k-e  Model  ....  1 6 1 

5.30  Turbulence  Intensities  for  Hobbs  Cascade  .  1 62 


5.31  Pressure  Distribution  for  Hobbs  Cascade,  Alg.  Model  164 

5.32  Overall  Performance  for  the  PSU  Compressor  Rotor  .  165 

5.33  Convergence  History  for  the  Rotor  Computation  ...  167 

5.34  Pressure  Distribution  at  R=.587  for  the  Rotor  .  .  .  168 

5.35  Pressure  Distribution  at  R=.670  for  the  Rotor  ...  169 

5.36  Pressure  Distribution  at  R=.750  for  the  Rotor  .  .  .  170 

5.37  Pressure  Distribution  at  R=.832  for  the  Rotor  ...  171 


xi 


5.38  Streamwise  Boundary  Layer  Profiles  on  the  Suction 

Side  for  the  Rotor .  172 

5.39  Streamwise  Boundary  Layer  Profiles  on  the  Pressure 

Side  for  the  Rotor .  174 

5.40  Radial  Velocity  Profiles  at  S=.99  for  the  Rotor  .  .  175 

5.41  Hub  Wall  Boundary  Layer  Profiles  for  the  Rotor  .  .  176 

5.42  Secondary  Flow  Vectors  at  X=.35  for  the  Rotor 

a)  Computed .  177 

b)  Measured:  Murthy  and  Lakshminarayana  (1987)  .  .  178 

5.43  Secondary  Flow  Vectors  at  X=.58  for  the  Rotor 

a)  Computed .  179 

b)  Measured:  Murthy  and  Lakshminarayana  (1987)  .  .  180 

5.44  Secondary  Flow  Vectors  at  X=1.0  for  the  Rotor 

a)  Computed .  1 8 1 

b)  Measured:  Murthy  and  Lakshminarayana  (1987)  .  .  182 

5.45  Pressure  Loss  Contours  at  the  Rotor  Exit 

a)  Computed .  1 84 

b)  Measured:  Murthy  and  Lakshminarayana  (1987)  .  .  185 

A . 1  Algebraic  3-D  Periodic  H-Grid 

a)  0-Z  plane  at  mid-span .  199 

b)  R-0  plane  at  mid-chord  .  200 

c)  R-Z  plane  at  mid-pitch .  201 


X 


LIST  OF  SYMBOLS 

A,B,C 

DfS  vectors  in  the  governing  equations 

An  Fourier  coefficient  matrix 

a  pressure  relaxation  coefficient 

c  chord  length 

p 

Cp  coefficient  of  pressure  2(p-pi)/oQi 

I  imaginary  unit  -fT\ 

J  Jacobian  of  the  grid  transformation 

k  turbulence  kinetic  energy 

ksec  mass  averaged  secondary  flow  kinetic  energy 
n  normal  direction 

p  static  pressure 

p0  stagnation  pressure  p  +  1/2pQ2 

Pk  production  of  turbulence  energy 

Q  total  velocity  vector 

q  dependent  vector  (p,u,v,w)T 

R  hub  to  tip  ratio 

Re  Reynolds  number 

S  streamwise  distance  non-dimens ionalized  by  c 

s  cascade  blade  spacing 

Ti  turbulence  intensity  /~k/1.5 

u,v,w  non-dimensional  velocity  components 

i  t 

u  ,v  ,w  time-averaged  turbulence  velocity  fluctuations 
u»  friction  velocity 

Us  streamwise  velocity  component 

Un  transverse  velocity  component 


xiii 


Wr  radial  velocity  component 

x-y-z  Cartesian  coordinates 
3  partial  derivative  w-r-t  subscript 

e  turbulence  dissipation 

tc  von  Karman  constant  (  =  .41) 

X  stagger  angle,  eigenvalue 

# 

X  complex  conjugate  of  X 

p,v  molecular,  kinematic  viscosity 

vt  non-dimensional  eddy  viscosity 

£2  non-dimensional  rotation  rate 

4>  flow  coefficient  Qm/Qt 

5-n-5  body  fitted  coordinates 
p  density 

Y  space  to  chord  ratio  s/c 

6  boundary  layer  thickness 

0  harmonic  frequency 

011  streamwise  momentum  thickness 

^loss  pressure-loss  coefficient  2(p0j-p0)/(p  Qt2) 
subscripts 
B  bulk 

<£  centerline 

e  edge  of  shear  layer 

I  inlet 

J  Jacobian 

i,J,k  indices  in  £,n,£  directions,  respectively 

m  mean 

n  normal 


p  first  point  from  boundary 

ps  pressure  side 

ss  suction  side 

s  streamwise 

t  tip 

superscripts 

as  assumed 


c  computed 

m  iteration  index  at  station  i 

n  global  iteration  level 
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CHAPTER  1 


INTRODUCTION 


1 . 1  The  Problem 

The  flow  in  turbomachinery  is  very  complex.  The  description  of  the 
inviscid  flow  is  generally  all  that  is  required  to  determine  the  gross 
properties  of  the  geometry,  such  as  blade  loading.  Some  important 
features  of  the  flow,  however,  are  due  to  the  viscous  nature  of  the 
actual  fluid  flow,  e.g.,  pressure  losses  and  inefficiencies.  Some  of 
the  viscous  flow  phenomena  found  in  a  typical  axial  flow  turbomachine 
can  be  seen  in  Figure  1.1.  These  include  the  three-dimensional  blade 
boundary  layer,  the  three-dimensional  wake,  and  overturning  of  the  flow 
in  the  end-wall  region.  If  the  geometry  in  question  is  a  rotor,  the 
flow  is  further  complicated  by  the  tip  leakage  flow  and  the  interaction 
of  the  leakage  jet  and  the  main  flow. 

For  three-dimensional  incompressible  flow,  three  momentum  transport 
equations  and  one  mass  conservation  equation,  collectively  known  as  the 
Navier-Stokes  equation,  govern  a  continuous  flow.  The  Navier-Stokes 
equation  can  be  simplified  to  give  algebraic  relations  for  the  dependent 
variables  of  static  pressure  and  flow  velocity  or  can  be  written  as  a 
single  scalar  Poisson  equation  for  a  potential  function.  Such  forms 
cannot  capture  the  above  mentioned  viscous  flow  phenomena.  The  full 
Navier-Stokes  equation  is  a  second  order,  non-linear  partial 
differential  equation  that  cannot  be  solved  exactly.  Thus,  researchers 
have  long  used  numerical  methods  to  solve  the  equation  on  discrete 
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points  in  the  domain  under  consideration. 

There  are  several  methods  for  solving  numerically  the  Navier-Stokes 
equation  in  full  or  simplified  form.  All  have  their  advantages;  some 
are  economical  yet  less  accurate  while  others  are  very  accurate  but 
require  large  amounts  of  computer  time.  A  description  of  some  of  the 
more  prominent  methods  for  solving  the  three-dimensional  Navier-Stokes 
equation  follows. 

1.2  Review  of  Related  Studies 

The  number  of  different  algorithms  designed  to  predict  viscous  flows 
is  vast.  However,  only  methods  developed  to  solve  the  three-dimensional 
Navier-Stokes  equation  are  of  interest  here.  The  types  of  methods  can 
be  catagorized  as  time-marching  methods,  space-marching  methods, 
parabolic-marching  methods  and  fully  elliptic  methods.  Time-marching 
methods  are  so  called  because  the  governing  equations  are  integrated  in 
the  time  coordinate.  Likewise,  space-marching  methods  use  a  streamwise 
integration  procedure.  Parabolic-marching  methods  are  essentially 
multiple  pass  versions  of  space-marching  methods  whereby  the  domain  is 
computed  several  times  to  relax  residual  errors.  Finally,  elliptic 
methods  are  relaxation  techniques  or  inversion  techniques  which  solve 
the  entire  matrix  generated  by  the  discretized  equations. 

1.2.1  Time-Marching  Methods 

The  unsteady  form  of  the  incompressible  Navier-Stokes  equation  seems 
well  suited  to  a  time  integration  procedure  except  that  there  is  no  time 
dependent  term  in  the  continuity  equation.  If  the  compressible  form  is 
used,  the  time- integration  can  proceed.  Pulliam  and  Steger  (1980)  used 


the  implicit  approximate  factori2ation  scheme  of  Beam  and  Warming  (1978) 
to  solve  the  compressible  equation  in  three  dimensions.  Even  though  the 
important  viscous  terms  were  retained,  the  method  required  fourth  order 
implicit  and  explicit  artificial  dissipation  for  convergence  to  be 
achieved.  The  method  was  used  to  compute  various  high  Reynolds  number 
supersonic  and  transonic  flows.  Computed  pressure  distributions 
compared  favorably  with  the  experimental  data.  If  one  were  to  use  the 
above  method  to  solve  very  low  Mach  number  flows,  difficulties  would  be 
encountered.  At  Mach  numbers  below  0.3  the  energy  equation  tends  to 
become  weakly  coupled  to  the  momentum  equations.  This  introduces  large 
numerical  errors  and  convergence  of  the  system  is  not  ensured.  Briley, 
Buggeln,  and  McDonald  (1985)  replaced  the  energy  equation  with  a 
constant  enthalpy  equation.  An  adiabatic  equation  of  state  was  used  to 
eliminate  the  pressure  as  a  dependent  variable.  The  resulting  system  is 
suitable  for  solving  low  Mach  number  flows.  The  Beam  and  Warming  (1978) 
algorithm  was  used  to  solve  the  coupled  system  for  a  laminar  horseshoe 
vortex  generated  by  a  wing-body  juction.  The  inlet  flow  Mach  number  was 
0.1.  Only  small  amounts  of  artificial  dissipation  were  required  to 
achieve  convergence.  Also  used  to  predict  three-dimensional 
compressible  flows  is  the  MacCorraick  implicit  method  (1982)  which  is 
based  on  the  author's  explicit  predictor-corrector  method.  It  only 
requires  the  solution  of  bidiagonal  matrices  and  so  is  very  efficient. 
Recently,  Dawes  (1986)  used  a  time-marching  method  to  solve  the  flow  in 
a  compressor  rotor.  A  finite  volume  formulation  was  used  to  discretize 
the  governing  equations.  Computed  Mach  number  contours  compared  well 
with  the  available  experimental  data. 

The  above  methods  use  the  compressible  formulation  which  assumes  a 


perfect  gas  is  being  considered.  If  the  fluid  is  water,  a  truly 
incompressible  form  of  the  equations  must  be  used.  Since  no  pressure 
term  is  present  in  the  continuity  equation,  solution  of  the 
incompressible  formulation  is  very  difficult  and  the  governing  equations 
must  be  manipulated  to  yield  a  more  tractible  form.  Following  the  ideas 
of  Chorin  (1967),  Kwak  et  al.  (19814)  introduced  a  time  derivative  of  the 
pressure  into  the  continuity  equation.  For  a  steady  flow,  a  time 
integration  of  the  system  yields  a  divergence  free  velocity  field.  The 
coupled  system  was  solved  using  the  Beam  and  Warming  algorithm  for 
several  three-dimensional  viscous  flow  problems  including  the  flow  in 
the  space  shuttle  main  engine  power  head.  Both  implicit  and  explicit 
fourth  order  artificial  dissipation  were  necessary  to  achieve 
convergence  which  was  seldom  greater  than  two  orders  of  magnitude. 

In  addition  to  the  above  implicit  methods,  explicit  methods  have 
been  used  to  compute  the  three-dimensional  viscous  flow.  Explicit 
schemes  are  generally  limited  to  very  small  time  steps  due  to  the  CFL 
condition  and  with  this  comes  large  computation  times.  However,  with 
the  advent  of  supercomputers,  explicit  methods  are  again  finding  favor 
among  researchers.  Shang  et  al.  (1980)  used  MacCormack's  explicit 
predictor-corrector  scheme  to  study  the  shock/boundary  layer  interaction 
in  a  wind  tunnel  diffuser.  The  explicit  scheme  allowed  for  a  high 
degree  of  vectorization  and  extremely  fast  computation  times  were 
reported.  Results  compared  reasonably  well  with  the  experimental  data. 
Chima  (1986)  used  a  vectorized  two-step  Runge-Kutta  scheme  to  solve 
quasi-three-dimensional  flows.  A  multigrid  algorithm  was  used  to 
efficiently  compute  the  flow  in  a  centrifugal  impeller.  Computed 
pressure  distributions  compared  well  with  the  experimental  data. 


6 


Artificial  dissipation  was  required  to  ensure  convergence.  Convergence 
to  five  orders  of  magnitude  was  enhanced  by  a  factor  of  three  when  using 
the  multigrid  scheme  over  a  single  grid. 

1.2.2  Space-Marching  Methods 

For  some  flows,  solving  the  full  elliptic  form  of  the  Navier-Stokes 
equation  is  not  necessary.  If  a  main  flow  direction  can  be  identified, 
the  governing  equation  can  be  "parabolized"  and  integrated  in  the  time¬ 
like  streamwise  direction.  The  Navier-Stokes  equation  is  parabolized  by 
neglecting  the  streamwise  diffusion  of  momentum  and  manipulating  the 
streamwise  flux  vector  such  that  the  product  of  the  inverse  of  its 
Jacobian  and  the  Jacobian  of  the  transverse  flux  vector  has  real 
eigenvalues.  In  supersonic  flow,  the  eigenvalues  are  all  real  and 
positive  if  the  streamwise  velocity  is  non-negative.  In  the  boundary 
layer,  one  eigenvalue  becomes  imaginary  and  the  system  becomes  unstable 
for  space-marching.  In  the  boundary  layer,  retaining  the  streamwise 
pressure  gradient  in  the  implicit  part  of  the  streamwise  flux  vector 
yields  an  unstable  system  as  the  grid  system  is  refined.  Such  departure 
solutions  are  well  documented.  Keeping  the  streamwise  pressure  gradient 
term  explicit  removes  this  instability  but  cannot  account  for  any  global 
interaction  between  the  viscous  and  inviscid  flows.  Shiff  and  Steger 
(1979),  following  boundary  layer  theory,  replaced  the  pressure  in  the 
viscous  layer  with  the  pressure  in  the  outer  supersonic  flow.  Vigneron, 
Rakich  and  Tannehill  (1978)  used  a  similar  marching  procedure,  namely 
Beam  and  Warming  (1978)  scheme,  yet  split  the  streamwise  pressure 
gradient  into  implicit  and  explicit  parts.  A  weighted  average  was  used 
based  on  the  flow  Mach  number.  Such  a  formulation  can  capture  mild 
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pressure  interaction  between  the  inner  and  outer  layers. 

Govindan  and  Lakshminarayana  (1987)  later  modified  Shiff  and 
Steger's  (1979)  method  for  internal  subsonic  flows.  A  special  bulk 
pressure  correction  was  required  to  conserve  global  mass  flow.  The 
method  proved  to  be  very  sensitive  to  the  value  of  the  explicit  assumed 
pressure  gradient  and  had  difficulty  in  low  Mach  number  flows. 

Along  very  different  lines,  Briley  and  McDonald  (1979)  based  a  new 
space-marching  method  on  secondary  flow  theory.  They  split  the  velocity 
into  four  distinct  parts  namely,  a  known  potential  velocity,  a  potential 
secondary  velocity,  a  solenoidal  secondary  velocity,  and  a  viscous 
streamwise  velocity  correction.  The  governing  equations  were  recast  to 
give  equations  for  each  new  component  of  the  velocity.  Elliptic  effects 
were  transmitted  through  the  potential  flow  which  was  required  a  priori. 
The  need  for  the  potential  flow  a  priori  is  a  major  drawback  for  the 
method  although  it  has  been  used  with  good  success  on  various  curved 
duct  geometries.  Along  similar  lines,  Dodge  (1976)  split  the  velocity 
into  potential  and  rotational  parts.  The  potential  component  was  used 
to  correct  the  pressure  so  that  continuity  was  satisfied. 

Patankar  and  Spalding  (1972)  developed  the  precusor  to  contemporary 
parabolic-marching  methods.  Their  space-marching  method  for  internal 
flows  could  solve  with  equal  success  the  compressible  or  incompressible 
Navier-Stokes  equation.  The  method  used  in  the  cross-plane  was 
essentially  the  SIMPLE  procedure  which  consisted  of  decoupling  the 
continuity  and  momentum  equations  and  solving  for  the  velocity 
components  on  a  staggered  grid  in  the  cross-plane.  The  streamwise 
momentum  equation  was  solved  for  the  streamwise  velocity  using  an 
assumed  streamwise  pressure  gradient.  The  transverse  momentum  equations 
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were  solved  in  turn.  A  bulk  passage  averaged  pressure  correction  was 
also  included  so  that  the  global  mass  flow  constraint  was  satisfied.  A 
two-dimensional  pressure  correction  was  also  included  in  order  to 
correct  the  velocities  to  satisfy  the  continuity  equation.  This  was 
accomplished  by  relating  the  gradients  of  the  pressure  correction  to  the 
velocity  corrections.  In  this  way,  a  two-dimensional  Poisson  equation 
for  the  pressure  correction  could  be  written  in  place  of  the  continuity 
equation.  Once  converged  at  a  particular  streamwise  station,  the 
pressure  correction  tended  to  zero  and  thus  the  velocity  corrections 
tended  to  zero  giving  a  divergence  free  velocity  field.  The  process  was 
iterated  to  convergence  at  each  streamwise  station  before  the  procedure 
was  advanced  to  the  next  station.  The  marching  scheme  could  be  used 
only  for  parabolic  flows  yet  was  a  major  breakthrough. 

Briley  (1974)  and  later  Ghia  and  Sokhey  (1977)  used  the  same 
pressure  split  scheme  and  similar  approaches  to  solve  the  incompressible 
equations  on  a  regular  grid.  Ghia  and  Sokhey  (1977)  introduced  an  extra 
inviscid  mean  pressure  which  was  a  function  of  the  streamwise  direction 
only.  The  inviscid  mean  pressure  was  assumed  known  and  its  presence 
allowed  for  the  computation  of  mildly  elliptic  flows.  Conservation  of 
global  mass  flow  was  used  to  determine  the  viscous  mean  pressure  and  the 
cross-plane  pressure  was  computed  from  a  Poisson  equation  which  was 
generated  by  taking  the  divergence  of  the  transverse  momentum  equations. 
Most  importantly,  the  system  did  not  in  general  satisfy  continuity. 
Thus,  irrotational  velocity  corrections  were  introduced  to  satisfy 
continuity.  Being  irrotational,  a  potential  function  could  be  written 
for  the  velocity  corrections  and  when  substituted  into  the  continuity 
equation,  yielded  a  second  Poisson  equation  to  be  solved  in  the  cross- 


plane.  Thus,  at  a  particular  streamwise  station,  the  streamwise 
momentum  equations  were  solved  for  their  respective  velocities  and  the 
pressure  was  corrected  to  satisfy  global  mass  flow  through  iterations. 
The  cross-plane  pressure  was  computed  from  the  Poisson  equation  and  the 
velocities  were  corrected  to  satisfy  continuity  by  solving  the  second 
Poisson  equation.  The  process  was  repeated  until  convergence  was 
achieved  and  the  procedure  was  advanced  to  the  next  streamwise  station. 
With  so  many  Poisson  equations,  the  method  proved  to  be  cumbersome  and 
thus  not  widely  used. 

To  remove  the  need  for  Poisson  equations,  Pouagare  and 
Lakshminarayana  (1986)  kept  the  pressure  gradients  implicit  and  solved 
the  entire  incompressible  formulation  in  a  coupled  fashion  using  the  LBI 
scheme  of  Briley  and  McDonald  (1980).  Thus,  local  continuity  and  the 
global  mass  flow  constraint  were  satisfied  without  the  need  for  solving 
extra  equations.  An  eigenvalue  analysis  showed  that  the  system  could  be 
space-marched  if  the  coefficient  of  the  streamwise  pressure  gradient  was 
less  than  zero.  This  result  could  be  inferred  from  the  work  of 
Vigneron,  Rakich  and  Tannehill  (1978)  for  a  Mach  number  of  zero.  Thus, 
the  streamwise  pressure  gradient  was  split  into  implicit  and  explicit 
parts  with  the  implicit  part  multiplied  by  a  small  egative  coefficient. 
For  small  values  of  the  coefficient,  the  method  was  relatively 
insensitive  to  the  explicit  assumed  pressure  gradient  when  applied  to 
parabolic  flows.  This  is  not  suprising  since  the  streamwise  pressure 
gradient  is  closely  tied  to  the  conservation  of  global  mass  flow.  The 
implicit  pressure  gradient  corrected  the  explicit  part  in  order  to 
conserve  global  mass.  The  method  was  very  fast  and  extremely  accurate 
for  straight  and  curved  ducts;  however,  several  problems  were 
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encountered  when  the  method  was  used  in  a  multi-pass  mode.  These 
problems  will  be  discused  extensively  in  Chapter  II. 

Recently,  due  to  the  speed  of  supercomputers,  explicit  methods  have 
been  used  to  solve  the  parabolized  Navier-Stokes  equation.  Spradley  and 
Stalnaker  (1981)  developed  an  interesting  procedure  which  used  the  time- 
dependent,  parabolized  Navier-Stokes  equation.  The  authors  termed  this 
a  Quasi-Parabolic  system  and  used  their  General  Interpolants  Method 
(GIM)  to  discretize  the  equation.  At  a  particular  streamwise  station, 
using  the  solution  at  the  previous  station  as  initial  conditions,  the 
equation  was  iterated  to  convergence  using  MacCormack's  predictor- 
corrector  scheme  in  the  cross-plane.  At  convergence,  the  solution 
procedure  was  advanced  to  the  next  streamwise  station  until  the  entire 
domain  was  computed.  Gielda  and  McRae  (1986)  solved  the  steady, 
parabolized  Navier-Stokes  equation  using  MacCormack's  two  step 
predictor-corrector  scheme,  marching  in  the  streamwise  direction.  Small 
marching  steps  were  required  due  to  CFL  conditions  yet  no  minimum  step 
size  requirement,  common  to  implicit  schemes,  had  been  encountei^d.  Due 
to  its  high  degree  of  vectorization,  the  method  proved  to  be 
computationally  efficient  even  with  the  fine  mesh  used.  It  should  be 
noted  that  the  fine  grid  resulted  in  increased  accuarcy  for  hypersonic 
flow  over  ogive  cylinders  at  an  angle  of  attack.  Also,  less  numerical 
damping  was  required  compared  to  implicit  methods. 

1.2.3  Parabolic-Marching  Methods 

Pratap  and  Spalding  (1976)  extended  Patankar  and  Spalding's  (1972) 


scheme  to  allow  for  a  multi-pass  solution  procedure.  The  domain  was 
computed  several  times  using  the  previously  computed  pressure  field  as  a 


new  assumed  pressure.  The  procedure  was  repeated  until  the  pressure 
corrections  dropped  to  zero.  In  this  way,  more  strongly  elliptic  flows 
could  be  computed  since  downstream  pressure  effects  could  be  transmitted 
upstream. 

Moore  and  Moore  (1979)  used  a  very  similar  procedure  to  compute  the 
viscous  compressible  flow  in  ducts.  A  finite  volume  method  was  used 
with  the  pressure  corrections  located  at  the  center  of  the  control 
volumes.  The  equations  were  integrated  in  the  streamwise  direction  and 
a  one-dimensional  pressure  correction  was  used  to  conserve  global  mass 
flow.  A  three-dimensional  pressure  correction  equation  was  written 
based  on  force  residuals  found  in  the  momentum  equations.  This  pressure 
correction  equation  was  solved  after  each  integration  of  the  momentum 
equations  until  the  force  residuals  dropped  to  zero.  Moore  and  Moore 
(1981)  later  improved  the  method  by  using  only  one  three-dimensional 
pressure  correction.  The  gradients  of  these  corrections  were 
manipulated  in  such  a  way  that  large  changes  in  the  transverse  gradients 
did  not  yield  large  changes  in  the  streamwise  gradient.  This  method 
could  be  used  on  non-orthogor.al  grids  whereas  the  earlier  method  had 
great  difficulty  with  such  grids.  The  three-dimensional  pressure 
correction  equation  based  on  force  residuals  was  retained  and  the  method 
was  used  to  compute  the  flow  in  a  ghost  impeller.  The  computed  static 
pressure  compared  well  with  the  experimental  data,  however,  the  computed 
velocity  field  did  not  compare  well. 

Rhie  (1983)  extended  Pratap  and  Spalding's  (1976)  multiple  pass 
marching  scheme  to  generalized  coordinates.  The  method  employed  a 
finite  volume  integration  and  pressure  corrections  to  conserve  local  and 
global  mass.  Averaging  was  used  on  control  volume  boundaries  to  give  a 
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finite  difference  expression  much  like  the  SIMPLE  scheme.  Numerical 
diffusion  was  introduced  to  avoid  instabilities  associated  with  central 


differencing  of  the  convective  terms.  The  computational  procedure  is 
very  similar  to  Pratap  and  Spalding's  procedure  except  that  once  a 
forward-marching  pass  was  completed,  a  three-dimensional  elliptic 
pressure  correction  equation  was  solved  to  accelerate  the  propagation  of 
downstream  pressure  corrections  upstream.  The  method  was  used  to 
compute  the  viscous  flow  in  curved  ducts  and  diffusers  and  in  impellers 
by  Rhie,  Delaney,  and  McKain  (1984).  In  all  cases,  the  results  compared 
favorably  with  the  experimental  data.  Khalil  and  Weber  (1984)  used  a 
procedure  very  similar  to  Rhie's  accept  that  special  attention  was  paid 
to  the  satisfaction  of  a  compatibility  relation  when  solving  the 
pressure  correction  equation.  The  compatibility  relation  is  based  on 
Green's  divergence  theorem  and  must  be  satisfied  exactly  if  a  converged 
solution  is  to  be  obtained  for  any  Poisson  equation.  The  method  was 
used  on  various  curved  duct  flows  with  very  good  results. 

1,2.4  Elliptic  Methods 

Elliptic  methods  solve  the  Navier-Stokes  equations  with  no 
parabolizing  assumptions  made.  Usually,  a  finite  volume  form  of  the 
uncoupled  equations  is  employed  which  is  relaxed  until  the  residuals 
drop  to  some  small  level.  Hah  (1984)  solved  the  uncoupled  equations  on 
a  staggered  grid.  The  staggered  grid  is  generally  difficult  to  use  on  a 
curvilinear  coordinate  system,  however,  Hah  overcomes  these  dificulties 
by  using  a  quadratic  upstream  interpolation  scheme  and  a  skew  upwinding 
scheme.  A  complex  algebraic  Reynolds  stress  model  was  used  to  model  the 
effects  of  the  turbulence.  This  method  was  used  to  compute  the  viscous 
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flow  in  a  turbine  end-wall  cascade  and  recently  in  a  compressor  rotor 
with  tip  clearance  (Hah  (1986)).  Results  compared  extremely  well  to  the 
experimental  data  for  all  test  cases.  Velocity  profiles  and  flow  angles 
compared  almost  exactly  to  the  experimental  data  in  the  tip  clearance 
region  even  with  as  few  as  five  grid  points  in  the  clearance  region. 
Moore  and  Moore  (1985)  used  nearly  the  same  procedure  as  Hah  but 
incorporated  an  algebraic  eddy  viscosity  model  for  the  turbulence.  The 
predicted  pressure  losses  were  in  better  agreement  with  the  experimental 
data  than  Hah's  prediction  for  a  turbine  end-wall  cascade. 

Vanka  (1985)  solved  the  governing  equations  in  a  coupled  fashion 
using  a  direct  solver.  Such  a  solver  requires  large  amounts  of  computer 
storage  but  can  be  very  efficient.  The  method  was  used  on  a  strongly 
curved  duct  with  very  good  results. 

Recently,  Rhie  (1986)  developed  a  full  elliptic  solver  which  employs 
the  Pressure  Implicit  Split  Operator  (PISO)  concept.  Several  levels  of 
pressure  corrections  are  used  to  correct  for  mass  flow  imbalances  and 
the  method  can  be  used  for  all  Mach  numbers  were  the  flow  remains  a 
continuum.  The  density  is  treated  implicitly  in  the  pressure  correction 
procedure.  A  predictor-corrector  type  algorithm  was  employed  and  a 
multigrid  procedure  was  included  to  enhance  convergence.  The  method  was 
tested  on  a  wide  range  of  flows  including  a  three-dimensional  driven 
cavity  and  a  turbine  end-wall  cascade.  Results  agree  very  well  with  the 
experimental  data. 


1.3  Objectives  and  Method  of  Approach 


The  major  objective  of  this  work  is  to  develop  a  new  parabolic¬ 
marching  method  to  solve  the  three-dimensional  incompressible  Navier- 
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Stokes  equation.  The  method  must  be  stable,  convergent,  and  accurate 
for  the  computation  of  a  wide  range  of  turbomachinery  flows.  The  new 
method  will  be  used  to  compute  blade  and  end  wall  boundary  layers  and 
three-dimensional  wakes  in  rotating  and  stationary  blade  rows.  A 
secondary  objective  is  to  discuss  drawbacks  of  a  few  earlier  methods  and 
to  include  high  order  turbulence  models.  Also,  some  of  the 
idiosyncracies  of  forward-marching  the  Navier-Stokes  equation  will  be 
discussed. 

In  order  to  develop  a  new  method,  the  difficulties  associated  with 
existing  techniques  must  be  discussed.  A  technique  recently  developed 
to  overcome  these  difficulties  will  be  analyzed  using  Fourier  stability 
theory  to  determine  its  global  stability  characteristics  as  a  multiple 
pass  method.  Next  a  simple  modification  to  the  technique  will  be 
analyzed  to  see  if  convergence  can  be  ensured  with  a  multiple  pass 
scheme.  Also,  an  eigenvalue  analysis  will  be  performed  to  determine  if 
the  modified  method  is  stable  during  the  forward-marching  integration 
process.  Finally  a  new  parabolic-marching  method  will  be  developed 
which  is  convergent  when  multiple  passes  are  made  and  which  overcomes 
some  of  the  difficulties  with  the  afore  mentioned  methods.  An 
eigenvalue  analysis  will  be  performed  to  ensure  stability  during  the 
forward-marching  process  and  a  Fourier  analysis  will  be  performed  to 
ensure  global  stability  during  the  iteration  process  for  a  wide  range  of 
flow  conditions.  A  new  solution  procedure  will  also  be  explored  which 
improves  predictions  for  coarse  grids.  The  new  method  employing  the  new 
solution  procedure  will  be  calibrated  by  computing  the  developing  flow 
in  a  straight  duct  with  square  cross-section.  Proper  transmission  of 
pressure  ellipticity  will  be  tested  by  computing  the  laminar  flow  in  an 
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S-shaped  duct.  Finally,  to  complete  the  validation,  various  turbulent 
turbomachinery  flows  will  be  computed  and  the  results  compared  to  the 
available  experimental  data. 
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CHAPTER  2 

ANALYSIS  AND  MODIFICATION  OF  EXISTING  TECHNIQUES 

2.1  Introduction 

The  equation  which  governs  fluid  flow  is  the  Navier-Stokes  equation. 
The  components  of  this  vector  equation  are  the  continuity  equation  and 
the  three  momentum  equations.  If  the  flow  is  compressible,  an  energy 
equation  is  added  to  the  set.  Only  incompressible  internal  flow  is 
considered  here.  Thus,  the  Navier-Stokes  equation  is  written  in  non- 
dimensional  form  as: 

+  dyB  ♦  —  1 /Re  (^xx^  ^  3yyD  +  "  S 

(2.1) 


where  Re  is  the  Reynolds  number  and  all  velocities  are  normalized  by 
some  bulk  inlet  velocity  and  the  pressure  by  the  dynamic  head. 

Equation  2.1  is  written  for  geometries  which  rotate  about  the  x  axis 
with  non-dimensional  rotation  R.  The  equation  is  second  order  elliptic 


due  to  the  diffusion  terms.  If  a  main  flow  direction  can  be  identified, 
as  with  most  turbomachinery  flows,  the  streamwise  diffusion  of  momentum 
can  be  neglected.  For  the  sake  of  example,  x  is  assumed  to  be  the 
streamwise  direction  and  so  neglect  the  diffusion  in  x  which  gives 
equation  2.2. 

3^A  +  3yB  +  3  ^ C  —  1/Re  ( 3yyD  +  9 z z ^ )  —  S 

(2.2) 

This  reduces  the  order  of  ellipticity  and  the  resulting  system  is 
termed  'partially  parabolic'.  With  such  a  system,  efficient  spacially 
marching  solution  procedures  may  be  used.  However,  the  streamwise 
pressure  gradient  term  introduces  first  order  ellipticity  and  thus 
requires  special  treatment. 

2,2  Difficulties  With  Existing  Techniques 

The  forward-marching  procedures  available,  whether  single  pass 
space-marching  or  multiple  pass  parabolic-marching  methods,  have  various 
drawbacks  associated  with  them.  The  following  methods  either  remove  the 
pressure  completely  from  the  formulation  or  split  the  pressure  into  a 
two-dimensional  transverse  pressure  and  a  one-dimensional  bulk 
streamwise  pressure.  These  modifications  are  useful  in  overcoming  the 
traditional  difficulties  associated  with  unmodified  pressure  terms. 
Briley  and  McDonald's  (1979)  method  removes  the  pressure  from  the 
formulation  and  requires  the  specification  of  the  potential  flow  for  the 
geometry  and  flow  conditions  under  consideration.  For  simple  straight 
and  constant  curvature  ducts,  the  potential  flow  is  almost  trivial. 
However,  for  complex  flows,  a  full  potential  solver  must  be  used 


initially  to  establish  the  potential  field.  This  may  require  a  large 
expenditure  of  computation  time  in  addition  to  the  viscous  flow  solver. 
Also,  the  viscous  pressure  field  is  not  a  computed  quantity  and  so 
viscous  pressure  losses  cannot  be  computed.  Patankar  and  Spalding's 
(1972)  method  solves  the  parabolized  Navier-Stokes  equation  on  a 
staggared  grid  in  an  uncoupled  fashion  using  a  pressure  split  scheme. 
The  staggared  grid  is  difficult  to  implement  in  complex  geometries  using 
generalized  coordinates.  Rhie  (1983)  later  recast  the  equations  for  a 
non-staggared  grid  but  like  Patankar  and  Spalding,  must  solve  a  Poisson 
equation  in  the  transverse  plane  for  the  pressure  corrections  in  order 
to  indirectly  satisfy  continuity.  Due  to  the  large  amount  of  computer 
time  required  to  reach  a  converged  solution,  a  Poisson  equation  should 
be  avoided  in  the  formulation.  For  this  reason,  Briley's  (1974)  and 
Ghia  and  Sokhey's  (1977)  methods  are  not  recommended,  in  fact  both 
methods  use  two  Poisson  equations  in  the  cross  plane,  one  for  the 
pressure  and  another  for  irrotational  velocity  corrections  which 
indirectly  satisfy  the  continuity  equation.  That  is,  the  continuity 
equation  in  the  above  methods  is  satisfied  only  through  an  iterative 
process  and  not  satisfied  directly.  Pouagare  and  Lakshminarayana  (1986) 
sought  to  remove  these  Poisson  equations  from  the  formulation.  To  do 
so,  the  flux  vectors  in  the  governing  equation  must  be  manipulated  to 
achieve  a  stable  forward-marching  procedure. 

In  the  following,  several  modifications  to  equation  2.2  are 
described  which  manipulate  the  pressure  terms  to  allow  for  stable 
forward-marching.  The  global  stability  of  multiple  passes  of  the  flow 
field  using  these  modifications  is  investigated  and  drawbacks  are 
described.  For  the  sake  of  example,  the  two-dimensional  form  of 
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equation  2.2  is  used  to  simplify  the  analysis,  however,  the  extension  to 
three  dimensions  is  straightforward  albeit  laborious. 

2.3  Method  of  Pouagare  and  Lakshminaravana  (PL) 

In  an  attempt  to  develop  an  efficient  single  pass  space-marching 
method  to  solve  the  partially  parabolized  Navier-Stokes  equation, 
Pouagare  and  Lakshminarayana  (1986)  modified  the  streamwise  pressure 
gradient  term.  The  modification  was  motivated  by  the  following 
eigenvalue  analysis  performed  by  them.  It  is  the  intent  of  this  author 
to  show,  through  a  Fourier  stability  analysis  in  section  2.3.2,  that 
this  method  is  globally  unstable  when  multiple  passes  of  the  flow  field 
are  made. 


2.3.1  Eigenvalue  Analysis 

For  simplicity,  the  following  analysis  is  performed  for  a  two- 
dimensional  system.  The  results  are  applicable  to  three  dimensions  as 
well.  First,  the  streamwise  pressure  gradient  in  the  flux  vector  A  is 
tagged  with  the  coefficient  a  giving: 


A= 


u 

2 

uc+op 

uv 


(2.3) 


The  two-dimensional  form  of  equation  2.2  is  then  linearized  and  written 
in  frozen  coefficient  form  as: 


Aj  3x9  ♦  3yq  -  Dj  3yy  q  =  S 


(2.4) 
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where  Aj,  Bj,  and  Dj  are  the  Jacobians  of  the  vectors  A,  B,  and  D, 

T 

respectively  and  q  is  the  dependent  vector  (p,  u,  v,  w)  . 
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1/Re 

0 

_  0 

0 

1/Re_ 

(2.-4  a) 


The  terms  on  the  left-hand  side  of  equation  2.4  and  all  other  equations 
are  treated  implicitly  in  the  computation  while  terms  on  the  right  are 
treated  explicitly.  For  stable  forward-marching,  the  eigenvalues  of  the 
matrix  Aj_1Dj  must  be  real  and  non-negative  for  proper  damping.  Also, 
the  eigenvalues  of  the  matrix  Aj-1Bj  must  be  real.  For  the  modified 
flux  vector  A  in  equation  2.3,  the  characteristic  equation  of  Aj"^Dj  is: 


\2(\  -  1/u  )  =  0 


(2.5) 


The  solution  of  equation  2.5  gives  the  following  eigenvalues: 


*1,2  =  0  >  *3  =  1/u 


(2.6) 


Therefore,  if  the  streamwise  velocity  is  positive,  the  scheme  is 
naturally  dissipative  and  the  first  condition  for  stable  forward- 
marching  is  satisfied.  For  Aj'^Bj,  the  characteristic  equation  is: 


\3  -  (v/u)X2  +  (1/o)X  -  v/(ou)  =  0 


(2.7) 


thus,  the  eigenvalues  are: 

X i  =  v/u  ,  X2,3  =  ±I//cT 


(2.8) 


where  I  =  /^T.  Therefore,  if  a  is  less  than  zero,  the  eigenvalues  are 
all  real.  In  such  a  case,  forward-marching  is  stable. 

Pouagare  and  Lakshminarayana  set  a  to  some  small  negative  value 
(-0.01).  With  this,  the  authors  modified  the  source  vector  S  for 
consistency  in  the  following  manner  (with  no  rotation): 


S  = 


0 

-3xpas  +  o3xp 
0 


(2.9) 


where  3xpas  is  an  assumed  pressure  gradient  and  is  forward  differenced 
for  proper  transmission  of  elliptic  effects  and  o3xp  is  differenced  in 
the  same  manner  as  the  implicit  part  namely: 


o3xp  =  o(pi  -  pi_i)/Ax 


(2. tO) 


where  i  is  the  streamwise  index  at  which  the  solution  is  desired.  Since 
S  is  a  source  term,  Pi  is  not  yet  known  and  so  must  be  replaced  by  pasi. 
Thus,  equation  2.10  becomes: 
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o3xp  =  o(paSi  -  pi_i)/Ax 


(2.11) 


The  final  system  is  inconsistent  unless  pas  =  p.  The  equation  was 
discretized  and  solved  using  the  LBI  scheme  of  Briley  and  McDonald 
(1980).  Pouagare  and  Lakshminarayana  used  a  very  small  o,  i.e.,  -0.01. 
This  relaxes  the  streamwise  pressure  gradient  condition  so  that  the 
global  mass  flow  constraint,  through  the  satisfaction  of  the  continuity 
equation  coupled  to  the  momentum  equations,  can  be  satisfied  directly. 


2.3.2  Global  Stability  Analysis 


With  the  modifications  described  in  section  2.3.1,  the  system  of 
equations  can  be  space-marched.  It  is  apparent  that  if  the  assumed 
pressure  field  is  not  correct,  the  computed  pressure  field  will  adjust 
so  that  the  governing  equations  will  be  satisfied.  One  must  now  ask 
whether  multiple  passes  of  the  domain  will  relax  the  residual  errors  due 
to  an  incorrect  assumed  pressure  such  that  p  -*■  pas.  The  following 
global  stability  analysis  will  answer  this  question. 

To  perform  the  analysis,  equation  2.2  must  be  discretized  in  the  way 
in  which  it  is  to  be  solved.  Two  point  backward  differences  are  used 
for  streamwise  derivatives  while  three  point  central  differences  are 
used  for  transverse  derivatives  giving: 

AJ  (qij-qi-lj)"/^  +  Bj  (qi,j+i-  Qi, j-1  )n/2Ay 
-  Dj  (qi,j+1-2qi,j+qi,j-l)n/(Ay)2  =  Sj  qn"1i,j 

(2.12) 


where  n  is  the  global  iteration  index  (pass  number)  and  J  is  the 


N 


transverse  index. 
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The  multiple  pass  stability  or  global  stability  of  the  method  can  be 
checked  using  a  von  Neumann  stability  analysis.  The  underlying  premise 
of  this  analysis  according  to  Smith  (1978)  is  that  the  error  function  at 
any  point  in  the  domain  satisfies  the  partial  differential  equation 
under  consideration.  Thus,  it  is  assumed  that  any  error  en£f j  satisfies 
equation  2.12.  The  error  is  also  separable  into  a  function  of  space  and 
a  function  of  iteration  level.  With  this,  the  error  can  be  decomposed 
in  a  Fourier  series.  The  necessary  and  sufficient  condition  for  the 
global  stability  of  a  two  level  system  is  that  the  error  due  to  the 
iteration  level  must  not  increase  with  increasing  iterations. 
Therefore,  the  following  relation  is  introduced  substituting  the  error 
function  for  the  dependent  vector  q  in  equation  2.12. 


9ni,J  -  eni,j  =  An  exP  i0x  ♦  J®y  ) 


(2.13) 


where  An  is  a  Fourier  coefficient  and  I  is  the  imaginary  unit.  Then, 
equation  2.12  becomes: 


Pj  An  =  Hj  An_ 1 


(2.14) 


where 


Pj 


0 

a(£i+£4) 

*2 


Hi 

2u£i+v£2~H3 

v£i 


*2 

U&2 

u£i+2v8-2"H3 


and 


24 


I 


H 


Hj 


0 

0+I-II5 

0 


0 

0 

0 


with 


ill  =  1  -  cos  9X  +  I  sin  9X 

£2  =  (Ax/Ay)  I  sin  0y 

£3  =  2Ax  (cos  0y  -  1 ) / { Re  (Ay)2} 

£.4  =  cos  0X  -  I  sin  0X 

£5  =  cos  8X  +  I  sin  0X 


The  amplification  matrix  G  is  now  defined  as: 


G  =  Pj-'  Hj  = 


Ml 
m2 

L  M3 


0  0 
0  0 
0  0 


where 


M1=(o+1-£5)(v£1£2-£i(u£1+2v£2-£3))/|Pj| 
M2=(o+1-£5)(-£22)/ j  Pj | 
M3=(o+1-£5)£1£2/|Pj| 
|Pj|=-£1(o(£1+£4)(u£1+2v£2-£3)-u£22)) 
+£2(o( £i+£4)v£i-£2(2u£i+v£2-£3) ) 


(2.15) 


(2.16) 


such  that 
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In  order  to  achieve  a  convergent  system,  the  absolute  value  of  the 
eigenvalues  of  G  must  be  less  than  or  equal  to  unity.  The  eigenvalues 
of  G  were  found  to  be: 

X-|  =  ( c+  I-II5)  {v8,iS,2-Hi  (uJl i+2 v8,2~®'3) }  /  I  I 

(2.18) 

X2,3  =  0 

(2.19) 

One  can  check  the  maximum  value  of  the  eigenvalues  by  setting  0x=9y=ir. 
This  gives: 


I  x Imax I  =  I  1  +  2/d  | 


(2.20) 


which  for  o=-0.01  (as  used  by  Pouagare  and  Lakshminarayana  (1986))  is 
199.  Note  that  equation  2.20  is  not  a  function  of  the  velocity  ratio 
nor  the  grid  aspect  ratio;  thus,  the  PL  method  is  unconditionally 
unstable  when  used  in  a  global  marching  procedure.  However,  it  is 
apparent  that  the  scheme  is  globally  stable  when  a  is  set  to  -1.  A 
close  inspection  of  the  PL  equations  reveals  that  this  gives  a  simple 
forward  difference  for  the  streamwise  pressure  gradient  which  is 
essentially  the  same  as  the  Rubin  and  Reddy  (1983)  scheme  on  a  regular 
grid. 

For  completeness,  the  global  stability  for  all  wave  numbers  must  be 
investigated.  Recall  that  if  the  absolute  value  of  X  must  be  less  than 
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« 

or  equal  to  unity  for  global  stability,  then  the  absolute  value  of  X  X 

ft 

must  also  be  less  than  or  equal  to  unity  where  X  denotes  the  complex 
conjugate  of  X.  With  u  and  v  set  to  one,  Ax  =  Ay  and  Re  arbitrarily 
set  to  1000,  the  global  stability  for  all  wave  numbers  for  o=-0.01  and 
-1.0  can  be  reviewed  in  Figures  2.1  and  2.2,  respectively.  The  hatched 
area  indicates  the  region  where  the  scheme  is  unstable.  One  can  see 
that  for  high  wave  numbers,  the  original  PL  method  with  o=-0.01  is 
unstable  but  with  o=-1.0,  the  scheme  is  stable  everywhere.  When  the 
velocity  ratio  is  changed,  there  is  no  change  in  the  stability 
characteristics.  There  is,  however,  a  large  change  when  the  grid  aspect 
ratio  is  changed.  The  effect  of  the  ratio  of  Ax/Ay  on  the  stability 
can  be  seen  in  Figures  2.3  and  2.4  for  o=- 0.01.  When  the  aspect  ration 
of  the  grid  is  changed  to  Ax/ Ay =0.01  from  unity,  the  system  of 
equations  becomes  ustable  for  a  wider  range  of  wave  numbers  but  the 
level  of  the  eigenvalues  remains  essentially  unchanged.  The  stable 
region  is  confined  to  wave  numbers  of  Gy/ir  <  0.04.  This  is  consistent 
with  the  well  known  stability  limitation  on  a  minimum  Ax  when  solving 
incompressible  flow  with  a  space-marching  method.  As  Ax  is  decreased, 
the  instability  becomes  greater.  Likewise,  when  Ax  is  increased, 
stability  is  enhanced  since  the  zone  of  ellipticity  is  overstepped. 
This  can  be  seen  for  Ax/Ay=100.  For  this  aspect  ratio,  there  appears 
to  be  a  wider  range  of  wave  numbers  where  the  global  iteration  procedure 
is  stable.  The  region  of  stability  is  confined  to  values  of  0x/tt 
between  0.1  and  0.9.  At  0x=0y=iT,  the  eigenvalue  still  reaches  its 
maximum  of  199.  If  the  absolute  value  of  the  eigenvalue  is  greater  than 
unity  for  any  wave  number  then  stability  cannot  be  assured.  Variation 
of  the  Reynolds  number  from  10  to  10^  produces  no  change  in  the  global 


Figure  2.1  Real  Part  of  X*X, 


Numbers  for  the  PL  method  ( 


Figure  2.4  Real  Part  of  X*X,  os-0.01  for  All  Have 
Numbers  for  the  PL  method  (u=v=1  ,Ax=100  Ay=1  ,Re=1000) 
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stability  characteristics.  Thus,  it  is  apparent  that  only  when  o=-1.0 
is  the  PL  scheme  stable  for  multiple  passes  of  the  flow  field. 
Unfortunately,  the  small  value  of  o  is  the  most  attractive  feature  of 
the  PL  method  since  this  allows  for  very  good  viscous  flow  predictions 
when  the  inviscid  pressure  is  given  as  the  assumed  pressure  (see 
Pouagare  and  Lakshminarayana  (1986)).  With  o=-1.0,  viscous  predictions 
of  complex  flows  become  difficult  since  the  streamwise  pressure  gradient 
condition  is  nolonger  relaxed  and  the  coupling  of  the  continuity 
equation  to  the  momentum  equations  stiffens  the  solution  process. 

2.3.3  Remarks  on  the  Method 

The  PL  method  was  found  to  be  relatively  insensitive  to  the  assumed 
pressure  distribution  for  strongly  parabolic  flows,  e.g.,  developing 
flow  in  a  constant  curvature  duct  with  constant  cross-sectional  area. 
This  is  no  surprise  since  the  satisfaction  of  the  global  mass  flow 
constraint  should  set  the  proper  streamwise  pressure  gradient  for  an 
internal  flow.  For  parabolic  flows,  the  transverse  pressure  gradients 
are  essentially  constant  or  change  only  slightly  in  the  streamwise 
direction.  Thus,  their  effect  can  be  introduced  as  some  source  term  in 
the  transverse  momentum  equations  leaving  the  pressure  as  a  dependent 
variable  in  only  one  equation,  the  streamwise  momentum  equation.  With 
continuity  satisfied  exactly,  the  proper  velocity  field  is  always 
computed  regardless  of  how  the  streamwise  pressure  gradient  is  treated. 
The  pressure  is  always  adjusted  during  the  computation  in  response  to 
the  assumed  pressure.  Combining  the  assumed  and  computed  pressure 
through  a  simple  functional  relationship  yields  an  integral  equation  for 
the  actual  pressure  as  was  done  by  Pouagare  and  Lakshminarayana  (1986). 
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When  the  transverse  pressure  gradients  vary  with  the  streamwise 
direction,  as  in  cascades,  the  actual  pressure  cannot  be  obtained  from  a 
simple  relationship.  When  using  the  PL  method  for  complex  flows,  while 
good  viscous  velocities  are  obtained,  the  viscous  pressure  cannot  be 
determined  easily.  One  alternative  is  to  use  the  Poisson  equation  for 
the  pressure  which  uses  only  the  viscous  velocity  field  to  compute  a  new 
viscous  pressure  field  (see  Kirtley  and  Lakshminarayana  (1985)).  The 
question  of  whether  the  Poisson  equation  can  be  coupled  to  the  PL  method 
to  yield  a  convergent  system  is  defered  to  section  2.5. 

The  method  is  very  useful  for  obtaining  the  viscous  velocity  field 
in  a  complex  geometry  if  a  fairly  accurate  presciption  of  the  inviscid 
pressure  field  is  used  as  the  assumed  pressure.  It  should  also  be  noted 
that,  due  to  the  small  value  of  a,  the  uncoupling  of  the  odd  and  even 
points  in  the  computation  is  severe. 

It  is  useful  to  comment  that  the  computation  of  complex  flows  is 
successful  with  the  PL  method  because  one  condition,  namely  the 
streamwise  pressure  gradient  condition,  is  relaxed  during  the 
computation.  It  should  also  be  noted  that  most  methods  described  in 
chapter  I  also  relax  some  condition,  usually  the  continuity  condition. 

Regarding  the  stability  of  the  multiple  pass  mode,  a  look  at  a  one¬ 
dimensional  system  as  suggested  by  Abdallah  (1987a)  leads  to  the  same 
conclusion  as  the  von  Neumann  analysis.  The  one-dimensional  system  is: 

3xu  =  0 

3xu2  +  o3xp  =  -3xpas  +  o3xp 

(2.21) 

where  the  second  term  on  the  right  side  of  the  momentum  equation  is 
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differenced  as  in  equation  2.11.  Inversion  of  the  system  shows  that: 
ui  =  ui^ 

Pi  =  l/o  Oxpas)  ♦  pasi 

(2.22) 

Thus,  regardless  of  the  value  of  pas,  the  correct  velocity  is  always 
computed.  On  the  otherhand,  the  difference  between  the  computed  and 
assumed  pressure  is  equal  to  the  assumed  pressure  gradient  divided  by 
o.  Thus,  any  error  in  the  assumed  pressure  gradient  will  be  multiplied 
by  a  very  large  number  since  a  is  very  small.  If  a  was  large,  the 
system  would  be  stable  but  physically  unrealistic.  Thus,  the  same 
conclusion  of  unconditional  instability  of  the  PL  method  can  be  made. 
One  wonders  whether  a  similar  modification  to  the  pressure  gradients 
could  be  made  in  such  a  way  as  to  give  results  as  good  as  the  PL  scheme 
in  single  pass  mode  yet  be  stable  in  multiple  pass  mode. 

2.4  Modified  PL  Method  (MPL) 

In  an  attempt  to  overcome  some  of  the  drawbacks  of  the  PL  technique, 
a  modified  PL  method  is  developed  in  this  section  based  on  the  pressure 
corrections  of  Moore  and  Moore  (1981).  The  pressure  is  split  into  an 
assumed  pressure  and  a  pressure  correction. 

P  =  Pas  +  Pc 

The  gradients  of  p  are  split  as  follows: 

3xP  =  *xP3S  +  3xP° 

9yP  =  3yPaS  +  b3ypC 


(2.23) 
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3zP  =  3zPaS  +  P3zP° 

(2.24) 

Where  as  Moore  and  Moore  solve  an  uncoupled  set  of  equations,  the 
present  MPL  method  solves  a  coupled  system  of  equations.  Thus,  in  two 
dimensions,  vectors  A,  B,  and  S  in  equation  2.4  become: 


u 

V 

0 

A= 

u2+pc 

B= 

UV 

S= 

-3xPaS 

.  uv 

_ v2+bpc. 

-  -3yPaS  . 

So  that  as  pc  -*•  0,  the  exact  equation  is  solved  within  the  accuracy  of 
the  discretization.  It  is  hoped  that  if  the  absolute  value  of  b  is 
large  then  the  global  instability  characteristic  of  the  PL  method  may 
not  arise  and  global  convergence  of  the  MPL  method  may  be  achieved. 
Moore  and  Moore  use  the  notation  1/e  which  is  equivalent  to  b  and  use  a 
value  of  e  equal  to  -0.01. 


2.4.1  Eigenvalue  Analysis 


The  stability  of  the  forward-marching  procedure  of  the  MPL  method 
can  be  studied  in  a  manner  similar  to  that  in  section  2.3.1.  The 
Jacobians  of  the  modified  vectors  A  and  B  are: 


0 

1 

0  ' 

'  0 

0 

1 

Aj= 

1 

2u 

0 

Bj= 

0 

V 

u 

.0 

V 

u  _ 

_  b 

0 

2v  _ 

Again,  for  stable  forward-marching,  the  eigenvalues  of  the  matrix  Aj-1Dj 
must  be  real  and  non-negative  for  proper  damping  and  the  eigenvalues  of 
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the  matrix  Aj-1Bj  must  also  be  real.  For  the  flux  vector  A  in  2.25,  the 
characteristic  equation  of  Aj-1Dj  is: 


X2(X  -  1/u  )  =  0 


(2.26) 


thus,  the  eigenvalues  are: 


*1,2  =  0  ,  X3  =  1/u 


(2.27) 


Therefore,  if  the  streamwise  velocity  is  positive,  the  scheme  is 
naturally  dissipative.  For  the  characteristic  equation  is: 


X^  -  (v/u)X2  +  bX  -  vb/u  =  0 


(2.28) 


Solving  equation  2.28  gives  the  following  eigenvalues: 
Xi  =  v/u  ,  X2f3  =  ±/^E 


(2.29) 


Therefore,  if  b  is  less  than  zero,  the  scheme  is  stable  for  forward¬ 
marching.  One  should  also  notice  that  the  eigenvalues  in  2.29  are  the 
same  as  those  of  the  PL  method  in  equation  2.8  if  b  is  equal  to  l/o. 
Therefore,  the  MPL  method  should  give  very  similar  results  to  those  of 
the  PL  method  if  b  is  set  to  -100.  So  now  there  is  a  large  coefficient 
multiplying  the  implicit  pressure  correction  gradient.  As  seen  from  the 
analysis  for  the  PL  method,  global  convergence  is  expected  with  the  MPL 
method  with  the  absolute  value  of  b  large. 
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2.4.2  Global  Stability  Analysis 


The  question  now  arises  whether  the  MPL  scheme  can  be  used  in  a 
multiple  pass  mode  to  reduce  residual  errors.  Using  the  procedure 
outlined  in  section  2.3.2  with  the  Fourier  decomposition  in  equation 


2.13  gives: 


Pj  An  =  Hj  An_1 


where 


(2.14) 


2u2,  -|+vfi,2-2.3 


uJ,1+2v2,2-8'3 


*2(b-1) 


This  gives  the  following  amplification  matrix  G: 


where 


(2.30) 


M1=(b-1 ) 42(ut1l2-i2(2ua1+v£2-a3) )/ | Pj | 
M2=(b-1)i1l22/|Pj| 


M3=(1-bH2*12/|Pj| 

|Pj|=-MMuli+2vl2-l3>-uba22)) 

+42(vH12-ba2(2uJL1+vS,2-)l3) ) 

and  !it2t3t4  are  defined  as  before  in  equation  2.15. 

The  eigenvalues  of  G  are  then: 

Xi  =  Mi 

*2,3  =  0 

(2.31) 

Previously,  the  maximum  value  of  X-|  was  found  for  9x  =  9y  =  ir.  For  this 
combination  of  wave  numbers,  X i  =0  which  would  indicate  that  the  system 

is  stable  for  multiple  passes.  This  is  not  the  case  as  a  plot  of  the 

* 

absolute  value  of  X  Xi  in  Figure  2.5  will  attest.  Again  the  hatched 
area  indicates  regions  of  instability.  For  a  value  of  b=-100,  u=v=1  and 
Re=1000,  the  values  of  the  eigenvalues  are  much  less  than  those  of  the 
PL  technique.  In  fact,  the  eigenvalues  are  only  slightly  higher  than 
unity  and  are  less  than  unity  for  most  wave  numbers.  However,  by 

n 

definition,  the  system  is  not  stable  if  |X  X |  >1  for  any  wave  number. 
As  with  the  PL  method,  the  stability  characteristics  remain  unchanged 
when  the  velocity  ratio  is  changed.  However,  when  the  grid  aspect  ratio 
is  changed  from  Ax/Ay=1  to  0.01,  it  seems  apparent  from  Figure  2.6  that 
the  scheme  is  globally  stable.  Unfortunately,  there  is  a  small  unstable 
region  confined  to  the  area  where  Qy/i\  =  0  at  which  location  the 
maximum  eigenvalue  equals  1.01  for  this  combination  of  parameters. 
Again,  this  indicates  that  global  stability  is  not  ensured.  As  with  the 
PL  method,  when  the  ratio  of  Ax/Ay=100  (see  Figure  2.7),  the  region 


stable 


Figure  2.5  Real  Part  of  X  X,  b=-100  for  All  Wave 
Numbers  for  the  MPL  method  (u=v=1  ,Ax=Ay=1  ,Re=1000) 


where  the  global  marching  procedure  is  stable  is  greater,  i.e.,  for  0.1 
<  0x/it  <  0.9.  Still,  the  maximum  eigenvalue  is  1.01  at  0x=ir  and  0y=O. 
Therefore,  the  MPL  method  is  unstable  for  a  multiple  pass  procedure.  It 
is  interesting  to  note  that,  for  the  same  combination  of  parameters,  the 
MPL  method  is  stable  for  the  combination  of  wave  numbers  where  the  PL 
method  is  unstable  and  vice  versa.  A  possible  explanation  for  this  is 
that  the  streamwise  momentum  equation  is  manipulated  in  the  PL  method 
whereas  the  transverse  momentum  equations  are  manipulated  for  the  MPL 
method . 

The  above  analysis  is  not  exact  since  the  actual  equation  is  not 
solved  in  frozen  coefficient  form  and  the  treatment  of  the  boundary 
conditions  has  not  been  included  in  the  analysis.  This  leads  one  to 
believe  that  the  system  may  be  stable  after  all.  Therefore,  it  seems 
worthwile  to  persue  this  technique. 

2.M.3  Remarks  on  the  Method 

The  MPL  method  appears  to  be  very  promising.  The  forward-marching 
eigenvalue  analysis  indicates  that  the  single  pass  results  may  be  as 
good  as  those  of  the  PL  technique.  Indeed,  it  can  be  used  efficiently 
as  a  single  pass  method  as  can  the  PL  method.  Unfortunately,  the  method 
is  not  stable  for  multiple  passes  of  the  domain  even  though  the  pressure 
gradient  is  multiplied  by  a  large  coefficent  which  seemed  to  be 
appropriate  for  global  convergence  as  mentioned  in  section  2.3.3.  A 
computation  of  the  developing  flow  in  a  straight  duct  with  a  square 
cross-section  shows  a  slow  divergence  when  multiple  passes  are  made. 
Figure  2.8  gives  the  convergence  history  for  this  test  case.  Due  to  the 
large  value  of  b,  the  same  uncoupling  of  odd  and  even  points  is  present 


as  with  the  PL  scheme.  This  also  leads  to  great  difficulty  in  computing 
spacially  periodic  flows.  The  scheme  has  been  used  for  a  zonal  equation 
method  developed  by  Warfield  and  Lakshminarayana  (1987).  In  their 
scheme,  the  MPL  method  is  used  in  parabolic  regions  of  the  flow  and  a 
time-marching  method  is  used  in  the  elliptic  regions  and  to  correct  the 
pressure  field  in  the  entire  domain.  With  another  equation  determining 
the  new  pressure  field,  the  MPL  method  in  essentially  a  multiple  pass 
mode  and  is  shown  to  be  convergent.  The  question  then  arises  whether 
both  the  PL  and  MPL  techniques  could  be  used  in  a  multiple  pass  mode  and 
be  stable  if  a  separate  equation  for  the  pressure  is  used  to  determine  a 
new  pressure  field  between  passes  of  the  domain.  The  obvious  choice  for 
the  separate  equation  is  the  Poisson  equation  for  the  pressure. 

2.5  Poisson  Equation  for  the  Pressure 

The  Poisson  equation  for  the  pressure  is  obtained  by  taking  the 
divergence  of  the  momentum  equation.  For  the  three-dimensional  system, 
the  equation  looks  like: 


*  P  =  h 


Ep  —  *(2(  3jjV3yU+3xw32U+3yw32v),*‘(  +(3yv)  +( ^2^^  ) 


(2.32) 


(2.33) 


There  is  no  analytical  solution  to  equation  2.32  since  the  source 
term  Ip  is  not  separable.  Equation  2.32  can  be  inverted  directly  for 
p,  however,  this  requires  immense  amounts  of  computer  storage  for  three- 
dimensional  computations.  Relaxation  techniques  are  the  most  efficient 
way  to  solve  2.32  but  a  convergent  solution  requires  the  satisfaction  of 
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a  compatibility  relation  which  arises  from  Green's  divergence  theorem, 


namely : 


/  /  dV  =  /  3np  dS 


(2.34) 


where  V  is  the  volume  of  the  domain,  S  is  the  area  enclosing  the  domain 
and  3n  is  the  derivative  normal  to  S. 

Equation  2.3 2  was  coded,  solved  using  a  semi-implicit  method,  and 
coupled  to  the  MPL  technique.  In  order  to  test  the  convergence  of  this 
system,  the  turbulent  flow  in  an  S-shaped  duct  was  computed.  The 
geometry  is  given  in  Figure  2.9  along  with  the  measurement  locations. 
The  flow  was  measured  by  Taylor  et  al.  (1982)  at  a  flow  Reynolds  number 
of  40000  based  on  the  duct  width.  More  details  about  this  particular 
flow  are  given  in  section  3.8.2. 

With  no  assumed  pressure  field,  the  MPL  method  could  not  compute 
past  the  second  bend  of  the  duct.  Here,  the  computed  flow  separated  and 
the  turbulence  model  introduced  instabilities.  Thus,  the  computation 
was  restarted  using  an  assumed  pressure  field  derived  from  a  sine 
function  with  the  measured  bulk  pressure  drop  imposed.  This  assumed 
pressure  was  both  smooth  and  close  to  the  experimental  data.  With  this, 
the  MPL  successfully  computed  the  entire  domain.  The  results  are 
presented  in  Figures  2.10  through  2.13. 

The  computed  streamwise  and  transverse  velocity  profiles  are 
compared  to  the  experimental  data  at  station  2  in  Figures  2.10  and  2.11, 
respectively.  Here,  the  streamwise  boundary  layer  thickness  is 
reasonably  well  predicted  but  the  secondary  flow,  while  having  the 
proper  sense,  is  not  large  enough.  This  may  be  due  to  inlet  conditions 
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MPL 

Data: 


Station  2,  MPL  method,  Re=40000. 


Station  2,  MPL  method,  Re=40000. 


Station  5,  MPL  method,  Re=40000 


Transverse  Velocity  Profiles  for  the  S-Duct 


with  a  lower  than  measured  vorticity  or  a  poor  resolution  of  the 
transverse  pressure  gradient.  The  streamwise  and  transverse  velocities 
are  again  compared  to  the  experimental  data  at  station  5  in  Figures  2.12 
and  2.13,  respectively.  Station  5  is  perceived  to  be  the  location  where 
elliptic  influences  are  the  greatest.  Near  the  side  walls,  the 
predictions  of  the  streamwise  velocity  are  not  good.  The  secondary  flow 
is  well  predicted  at  all  transverse  locations  except  at  Y/D=0.9.  At 
this  location,  almost  no  secondary  flow  is  computed  which  is  not 
consistent  with  the  data.  These  below  average  predictions  should  not  be 
worrysome,  however,  since  this  is  the  first  pass  of  what  is  hoped  to  be 
a  multiple  pass  procedure  in  which  the  solution  converges  to  the  correct 
solution  after  several  global  iterations  with  the  Poisson  equation  used 
to  update  the  pressure. 

The  viscous  pressure  was  computed  from  equation  2.32  and  is  given  in 
Figure  2. 14.  Neumann  boundary  conditions  based  on  the  momentum 
equations  were  used  on  all  boundaries  save  the  downstream  boundary  where 
a  Dirichlet  condition  was  used.  One  can  see  from  Figure  2.14  that  the 
new  computed  pressure  field,  while  qualitatively  correct,  does  not  match 
the  experimental  data  near  Y/D=0.1  and  0.5.  When  this  pressure  was  used 
as  the  new  assumed  pressure  and  the  marching  process  repeated  with  the 
MPL  technique,  the  computation  exhibited  divergent  behavior  near  the 
begining  of  the  second  bend.  Manipulation  of  the  boundary  conditions 
for  the  pressure  solver  did  not  improve  the  convergence.  When  the 
laminar  flow  was  computed,  which  has  much  stronger  secondary  flows  and 
pressure  gradients,  the  same  divergent  behavior  was  observed.  From 
this,  one  may  conclude  that  the  Poisson  equation  coupled  to  the  MPL 
method  as  outlined  above,  does  not  yield  a  globally  stable  system  of 


equations.  There  appears  to  be  no  mechanism  in  the  system,  of  equations 
to  relax  the  pressure  errors.  Another  problem  may  be  due  to  the 
following.  When  using  Neumann  boundary  conditions,  equation  2.34  is 
seldom  satisfied  by  the  discretized  form  of  equation  2.32.  This  leads 
to  a  divergent  solution  of  the  Poisson  equation.  If  one  Dirichlet 
boundary  condition  is  used  then  equation  2.34  will  be  satisfied  exactly 
and  a  convergent  solution  will  result.  This  is  because  3np  on  the 
boundary  with  the  Dirichlet  condition  will  adjust  such  that  equation 
2.34  is  satisfied  as  it  must.  Unfortunately,  this  pressure  gradient  may 
not  satisfy  the  normal  momentum  equation  with  the  velocity  field  used  in 
equation  2.33.  This  inconsistency  will  not  affect  the  relaxation 
solution  of  equation  2.32  but  may  infect  the  pressure  field  in  such  a 
way  that  a  global  iteration  procedure  between  the  PL  or  MPL  equation 
system  and  equation  2.32  will  not  converge. 

In  other  words,  3np  based  on  the  momentum  equations  cannot  be  used 
on  all  boundaries  because  the  compatibility  relation  is  not  satisfied  by 
the  discretized  form  of  equation  2.32.  A  Dirichlet  condition  on  one 
boundary  cannot  be  used  because  the  3np  arising  from  the  solution  of 
equation  2.32  and  satisfaction  of  the  compatibility  relation  will  not  be 
consistent  with  the  3np  as  determined  from  the  normal  momentum  equation 
with  the  computed  velocity  field. 

Recently,  Abdallah  (1987)  developed  a  discretization  scheme  for  the 
Poisson  equation  such  that  the  compatibility  condition  is  satisfied 
exactly  by  the  discretized  equations.  Equation  2.32  was  shown  to 
converge  using  Neumann  boundary  conditions  based  on  the  momentum 
equations.  The  Navier-Stokes  equation  used  in  this  analysis  was 
discretized  using  a  central  difference  scheme.  The  me. rod  has  not  been 
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extended  to  the  backward  difference  scheme  adopted  here  for  the  forward¬ 
marching  solution  of  the  Navier-Stokes  equation. 

2.6  Conclusion 

The  PL  method  was  shown  to  be  unconditionally  unstable  in  a  multiple 
pass  marching  procedure  with  a  equal  to  any  value  other  than  -1.0. 
With  os-1,  a  standard  forward  difference  scheme  for  the  streamwise 
pressure  gradient  results  which  has  been  shown  to  be  globally  stable  as 
other  investigators  have  found  (see  Rubin  and  Reddy,  1983).  The 
advantage  of  the  standard  PL  method  (o=-0.01)  is  that  complex  flows  can 
be  computed  with  a  single  pass  of  the  domain  provided  a  good  assumed 
pressure  is  given  a  priori.  The  success  can  be  directly  attributed  to 
the  fact  that  one  condition,  namely  the  streamwise  pressure  gradient 
condition,  was  relaxed,  due  to  the  small  value  of  a,  in  favor  of 
satisfying  the  global  mass  flow  constraint.  This  was  accomplished  by 
coupling  the  continuity  equation  directly  to  the  momentum  equations 
during  the  solution  process. 

The  global  stability  of  the  PL  method  suffered  due  to  the  small 
value  of  o,  thus,  a  modified  version  of  the  method  (MPL)  was  developed 
in  an  attempt  to  overcome  this  lack  of  global  stability.  Even  with  a 
large  coefficient  multiplying  the  pressure  gradients,  the  MPL  method  was 
found  to  be  globally  unstable.  When  the  Poisson  equation  for  the 
pressure  was  introduced,  global  stability  was  not  assured  as  the 
computations  proved.  It  should  be  noted  that  the  solution  of  the 
Poisson  equation  in  three  dimensions  requires  large  amounts  of  computer 
storage  and  computation  time.  This  seems  to  reduce  the  advantages 
gained  by  using  parabolic-marching  algorithms  over  time-marching 


54 

algorithms.  It  is  evident  that  a  new  parabolic-marching  method  is  in 
order  which  overcomes  the  most  important  hurdle,  that  of  global 
stability. 
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CHAPTER  3 


n 


NEW  PARABOLIC-MARCHING  METHOD  (NPM) 


3.1  Introduction 

It  was  found  in  the  preceeding  analysis  that  the  PL  method  was 
unconditionally  unstable  for  a  global,  multiple  pass,  marching  procedure 
for  all  o  other  than  -1.  However,  if  o  is  equal  to  -1,  a  simple 
forward  difference  for  the  streamwise  pressure  gradient  is  recovered 
with  the  downstream  pressure  known  from  a  previous  pass  of  the  domain. 
The  maximum  eigenvalue  of  the  amplification  matrix  is  unity;  thus,  the 
system  is  stable.  This  yields  a  method  very  similar  to  that  of  Rubin 

■  and  Reddy  (1983)  except  that  it  is  on  a  regular  grid.  The  problem  with 
such  a  method  is  that  it  is  very  difficult  to  compute  the  first  pass  of 
a  complex  flow  unless  a  very  good  prescription  of  the  pressure  field  is 

■  given  initially.  The  PL  and  MPL  methods  generally  have  little  trouble 
computing  the  first  pass  since  one  condition  is  relaxed  during  the 
marching  process.  Recall  that  in  the  PL  method  the  streamwise  pressure 
gradient  is  relaxed  and  in  the  MPL  method  the  transverse  pressure 
gradients  are  relaxed.  This  allows  for  a  more  stable  forward-marching 
process  with  the  continuity  equation  coupled  to  the  momentum  equations. 

What  is  apparent  is  that  the  pressure  gradients  cannot  be  relaxed  in 
order  to  achieve  a  multiple  pass  procedure  that  is  stable.  Yet  some 
condition  must  be  relaxed  during  the  computation  in  order  to  solve 
complex  flows.  The  only  condition  left  to  be  relaxed  in  a  practical  way 
is  continuity. 


The  method  of  Chorin  (1967),  developed  for  incompressible  flow, 
relaxes  the  continuity  constraint  while  keeping  the  equation  coupled  to 


the  momentum  equations.  In  Chorin 's  method,  an  artificial  time 
derivative  of  the  pressure  is  added  to  the  continuity  equation  as  in  the 
following: 


1/8  3tP  +  V*Q  =  0 


(3.1) 


where  Q  is  the  total  velocity  vector.  This  relaxes  the  continuity 
constraint  during  a  time-marching  integration  of  the  Navier-Stokes 
equation  (see  Kwak  et.  al  19811).  At  convergence,  for  steady  flows,  the 
computed  velocity  field  is  divergence  free  since  9^p  0. 

Since  a  parabolic-marching  scheme  is  under  consideration  here,  and 
no  time  derivatives  are  present,  the  continuity  equation  must  be  written 
in  a  form  other  than  that  in  equation  3.1.  Successive  passes  of  the 
domain  may  be  thought  of  as  an  advancement  in  time.  With  this  in  mind, 
the  continuity  equation  is  rewritten  as: 


a(pn-pn_1)  +  VQ  =  0 


(3.2) 


where  n  is  a  global  iteration  index  and  a  is  a  coefficient  analagous  to 
1/8At  in  Chorin' s  method  where  B  is  some  positive  constant.  As 
determined  from  the  eigenvalue  analysis  performed  in  section  2.3.1,  the 
streamwise  pressure  gradient  must  be  forward  differenced  for  stable 
forward-marching.  Therefore,  forward  differencing  of  the  streamwise 
pressure  gradient  and  the  replacement  of  the  continuity  equation  with 
equation  3.2  yields  the  following  form  of  the  governing  equation: 
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3XA  +  dyB  +  —  1/Re  ( 3yyD  +  ^zz^)  ~  S 

(3.3) 
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It  should  be  noted  that  without  the  implicit  pressure  term  in  the 
streamwise  flux  vector,  a  globally  unstable  system  results. 

3.2  Eigenvalue  Analysis 

The  question  now  arises  whether  equation  3.3  can  be  forward-marched 
in  a  stable  fashion.  Using  the  analysis  described  in  section  2.3.1,  for 
a  model  two-dimensional  problem,  the  flux  vectors  A  and  B  are  written 
as: 


ap+u 

v 

A= 

u2-P 

B= 

uv 

uv 

2 

L  v^+p 

(3.4) 


The  negative  sign  of  the  pressure  term  in  A  is  due  to  the  forward 
differencing  procedure  used  for  the  streamwise  pressure  gradient.  With 
this,  the  two-dimensional  form  of  equation  3.3  is  linearized  and  written 


in  frozen  coefficient  form  as: 


Aj  SxQ  +  3yyQ  =  S 


(3.5) 


where  Aj,  Bj,  and  Dj  are  the  Jacobians  of  the  vectors  A,  B,  and  D, 

T 

respectively  and  q  is  the  dependent  vector  (p,  u,  v,  w)  . 
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For  stable  forward-marching,  the  eigenvalues  of  the  matrix  Aj_1Dj  must 
be  real  and  non-negative  for  proper  damping.  In  addition,  the 

eigenvalues  of  the  matrix  Aj_1Bj  must  be  real.  For  the  modified  flux 
vector  A  in  equation  3.^,  the  characteristic  equation  of  Aj~^Dj  is: 

X3-X2{ 1/u+a/( 1+2ua) }+X( 1/u)a/( 1+2ua)=0 

(3.6) 


The  solution  of  equation  3.6  gives  the  folowing  eigenvalues: 
X  i  =  0  ,  X2  =  1/u  ,  X3=a/(1+2ua) 


(3.7) 


It  is  apparent  that  a  can  be  negative  and  still  give  a  system  which  is 
naturally  dissipative  if  u  is  positive.  Since  the  analogy  was  drawn  to 
time-marching  methods  where  a  looks  like  1/BAt,  a  negative  value  for  a 
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is  not  physically  realistic.  The  global  stability  analysis  performed  in 
section  3.3  indicates  that  a  must  be  non-negative  for  global  stability 
and  convergence.  For  Aj-1Bj,  the  characteristic  equation  is: 

X^(1+2ua)  -  X2{2av+v/u{ 1/( 1+2ua) )}  +  X(2av2/u-1 )+v/u  =  0 

(3.8) 


Solution  of  equation  3.8  gives  the  following  eigenvalues: 


X]  =  v/u 

X2,3  =  {av±/(av)2+(  1+2ua)2} /( 1+2ua) 


(3.9) 


One  can  see  that  the  radicand  of  equation  3.9  is  always  positive 
regardless  of  the  value  of  a;  thus,  the  eigenvalues  are  all  real  and 
forward-marching  will  be  a  stable  solution  technique. 

3.3  Global  Stability  Analysis 


The  motivation  for  making  the  above  modification  to  the  governing 
equation  is  to  develop  a  marching  scheme  which  is  stable  in  multiple 
pass  mode;  thus,  the  following  stability  analysis  must  be  performed. 
Using  the  Fourier  decomposition  in  equation  2.1 3  for  the  discretized 
equation  in  2.12  with  the  modifications  made  in  equation  3.4,  the 
following  relationships  result: 


Pj  An  =  Hj  An_1 


(2.14) 


A 


« 


where 
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and 


Pj  = 


a(  £1+8,4)  8,1  8,2 

-(8,i  +  8,4>  2u8,i+v8,2~8,3  u8,2 

j_  8,2  v8,i  u8,i+2v8,2-8,3 
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a 
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0 
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0 
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+  I 
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The  amplification  matrix  G  is: 

0 
0 
0 


G  =  Pj'1  Hj  = 


Ml 
M2 
L  m3 


0 

0 

0 


where 

M 1  =  ( a(  ( 2u  8, 1  +  v  £2- )  ( u  8, 1 +2  v  8,2  -  8,3 ) -uv  8, 1 8-2 ) 
+£5(£l(u£i+2v£2-£3)-v£i£2))/|Pj| 


(3.10) 


(3.11) 
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M2=(a(  ( 8,1  +  8.14)  (uii+2v£2- £3  )+u82^ ) 

-t5(a(  Hi+ll4)(u1l1+2vil2-Jl3)-Jl22) )/  |  Pj  | 

M3=(a(-(  Jl  1+8,4) v£i-2.2(2ull  1+VH2-S.3)  )  +  fi-5(a(  H  1+4.4 )v£  1-8,18,2) )/  |  Pj  | 

|  Pj  I  =a(ai  +  8,4)((2u4i+v42-a3)(u8,i+2v8,2-2'3)-uv2'l2’2) 

+2- 1  ( ( a  1 -*■  2-4 )  ( u8 1 +2v  82_2’3 ) +u  2>22 ) 

-a2( (ai+a4)vai+a2(2uai+2va2-a3)) 

The  eigenvalues  of  G  are  then: 

Xi  =  Mi 

X2)3  =  0 

(3.12) 


The  maximum  value  of  Xi  can  be  determined  by  setting  0x=0y=n  which 
gives: 


Xi  |  =  |  1-2/[a(2u+c)  +  1] | 


(3.13) 


where  c=2Ax/Re/( Ay)^ 

Therefore,  if  a  is  non-negative,  the  maximum  eigenvalue  is  always 
less  than  unity  regardless  of  the  value  of  c.  Thus,  the  scheme  appears 
to  be  unconditionally  stable  in  multiple  pass  mode.  From  equation  3.13, 
one  can  determine  that  the  lower  bound  for  a  is  zero.  For  a  less  than 
zero,  the  maximum  eigenvalue  is  greater  than  unity  and  the  error  will 
grow  without  bound.  This  is  not  surprising  since  a  is  analagous  to 
1/BAt  in  time-marching  methods.  For  completeness,  the  stability 
characteristics  for  all  wave  numbers  for  various  combinations  of 

parameters  should  be  studied.  Again  let  u=v=1  and  Re=  1000  and  plot  the 

* 

real  part  of  X  Xi  for  all  wave  numbers  in  Figure  3.1  for  a  arbitrarily 
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set  to  0.2.  With  this  set  of  parameters,  one  can  see  that  the  method  is 

|| 

stable  for  all  wave  numbers.  The  maximum  values  of  X  Xi  occur  at  the 
limits  of  the  frequency  spectrum.  The  minimum  values  are  found  near 
0y/x“O.6.  With  the  ratio  u/v  changed  from  unity,  the  stability 
characteristics  do  not  change.  Also  not  an  important  parameter  is  the 
Reynolds  number  which  was  changed  from  10  to  10^  with  no  significant 
difference  in  the  stability  characteristics.  With  the  ratio 
Ax/Ay=0.01,  the  stability  characteristics  are  plotted  in  Figure  3.2. 

It  is  evident  that  the  maximums  are  now  confined  to  0y/iT=O  and  1, 
however,  the  scheme  is  stable  everywhere.  The  minimums  are  still  found 

near  0y/x=:O.6.  With  Ax/Ay  set  to  100,  the  stability  characteristics 

* 

are  much  different  (see  Figure  3.3).  The  maximum  values  of  X  X]  are 
now  confined  to  a  very  thin  region  for  0x/tt  less  than  0.01  and  greater 
than  0.99.  Large  values  of  a  should  severely  underrelax  the  dependent 
variables  as  does  a  small  At  in  a  time-marching  method.  The  value  of  a 
was  varied  from  1  to  100  and  the  differences  in  the  global  stability 
characteristics  are  given  in  Figure  3.4.  One  can  see  that  there  is  no 
value  of  1  Xj  greater  than  unity  so  that  the  system  is  stable  for  a 
multiple  pass  solution  procedure. 

The  above  analysis  was  performed  for  a  two-dimensional  system  and  it 
was  assumed  that  the  results  were  applicable  to  three-dimensions.  For 
completeness,  the  stability  analysis  should  be  performed  for  three- 
dimensions.  Recall  that: 


Pj  An  =  Hj  An“ 1 


(2.14) 


where  Pj  and  Hj  are  now  given  as: 


.996 

.995 

.994 
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where  8-1,2  3,^,5  are  defined  in  equation  2.15  and 

8-6  =  (Ax/Az)  I  sin  0Z 

8-7  =  2Ax  (cos  02  -  1)/(Re(A2)2) 

For  0x  =  0y  =  0z  =  Tr  the  maximum  value  of  X-|  is 

|  X *1  |  =  |  1-2/[a(2u+c+d)+1  ]  | 

(3.13a) 

where  c=2Ax/Re/( Ay)2  and  d=2Ax/Re/(Az)2.  Thus,  if  a  is  positive,  the 
system  in  three-dimensions  is  stable  for  multiple  passes.  Since  the 
discretization  of  the  governing  equation  is  consistent,  by  Lax's  theorem 
one  may  state  that  the  system  is  also  convergent.  In  actual  practice, 
no  method  is  always  convergent  for  complex  turbulent  flows,  however, 
confidence  is  high  that  the  new  method  will  outperform  both  the  PL  and 


MPL  methods. 
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3.4  Convergence  of  the  Method 


It  is  important  to  show  that  the  method  is  able  to  converge  to 
machine  precision  to  show  that  artificial  energy  losses  are  not 
introduced  with  the  introduction  of  the  pressure  term  in  the  continuity 
equation.  The  developing  laminar  flow  in  a  straight  duct  with  square 
cross  section  was  computed  on  a  coarse  grid  with  a=0.01.  The  resulting 
convergence  history  is  given  in  Figure  3.5.  Convergence  to  machine 
precision  (1  x  10“^  on  the  IBM  3090-200)  is  achieved  and  no  divergent 
behavior  is  apparent.  Regarding  convergence  rates,  it  is  known  that  the 
value  of  the  maximum  eigenvalue  sets  the  convergence  rate  of  the  global 
iteration  process.  That  is, 


en/en-1  a  1/ | Xt  |  =  1/p(G) 


(3.14) 


where  p(G)  is  the  spectral  radius  of  the  amplification  matrix  G.  Thus, 
if  X -j  is  minimized,  the  rate  of  convergence  will  be  maximized.  It  is 
apparent  that  equation  3.13  can  be  minimized  if: 


a=  1/(2u+c) 


(3.14a) 


For  u=v= 1 ,  Re=  1000,  Ax=Ay=1,  the  optimum  value  of  a  is  0.5.  This  may 
reduce  the  largest  value  of  the  eigenvalue  but  comparison  of  Figures  3.1 
and  3-6  shows  that  for  other  wave  numbers,  the  eigenvalues  for  a=0.5  are 
generally  higher  than  for  a=0.2.  Therefore,  the  apparent  optimum  value 
of  a  may  not  yield  the  fastest  convergence.  Similarly  for  three- 
dimensions,  the  optimum  value  for  a  is: 


Figure  3.5  Convergence  His 
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a=  1/(2u+c+d) 


(3.14b) 


A  problem  with  parabolic-marching  methods  in  general  is  that,  since 
pressure  information  propagates  upstream  only  one  streamwise  step  per 
global  pass,  the  downstream  pressure  boundary  condition  will  not 
influence  the  inlet  flow  until  the  same  number  of  passes  as  streamwise 
steps  is  made.  Also,  small  wave  length  errors,  order  Ax,  are  damped 
out  quickly  while  large  wave  length  errors,  order  NAx,  are  damped  out 
at  a  much  slower  rate  rate  or  not  at  all.  Because  of  this,  convergence 
of  any  parabolic-marching  method  is  very  slow  and  the  NPM  method  is  no 
different . 

One  technique  for  accelerating  the  convergence  of  the  numerical 
method  is  the  multigrid  procedure  of  Brandt  (1977).  In  incorporating 
this  into  the  NPM  technique,  the  multigrid  procedure  would  be  applied  in 
the  streamwise  direction  only.  The  coarse  grid  equations  would  be 
solved  on  grids  with  a  progressively  increasing  aspect  ratio.  While 
this  does  not  reduce  the  stability  of  the  NPM  method  as  seen  in  section 
3.3,  it  may  severely  decrease  accuracy  if  the  prolongation  operation  of 
the  multigrid  is  not  performed  perfectly.  Rubin  and  Reddy  (1983)  have 
used  multigrid  with  their  parabolic-marching  method  but  the  improvement 
in  convergence  does  not  seem  to  be  worth  the  extra  computational  effort. 

Another  method  for  accelerating  the  convergence  of  a  parabolic¬ 
marching  method  was  recently  developed  by  Tenpas  and  Pletcher  (1987). 
Their  method  solves  a  modified  Poisson  equation  for  a  pressure 
correction  with  a  single  backward  pass  of  the  domain  after  each  forward¬ 
marching  pass  of  the  mean  flow  equations.  A  variable  coefficient 
multiplies  the  transverse  pressure  gradient  in  the  Poisson  equation 
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which  is  different  for  each  case.  The  increase  in  convergence  is 
impressive  but  results  for  a  flow  with  large  transverse  pressure 
gradients  have  not  been  published  by  them  and  one  must  wonder  what  the 
influence  is  of  the  modified  transverse  pressure  gradient  in  the  Poisson 
equation. 


3.5  Transformation  and  Discretization  of  the  Equation 


In  order  to  implement  boundary  conditions  in  a  straightforward 
fashion  for  complex  geometries,  it  is  desirable  to  solve  the  governing 
equation  on  a  generalized  coordinate  system.  Therefore,  the  equations 
solved  must  be  recast  into  the  generalized  coordinates  using  the 
following  transformation: 


5  =  5(x,y,z) 
n  =  n(x,y,z) 
c  =  £(x,y ,z) 

The  full  equation,  2.1,  is  transformed  using  a  chain  rule 
conservative  formulation.  The  Cartesian  velocity  components  and  the 
pressure  are  still  the  dependent  variables,  however,  the  equation  is  now 
solved  on  a  generalized  coordinate  system.  The  modifications  made  to 
equation  2.1  in  the  x-direction  which  yielded  equation  3.3  are  now  made 
in  the  streamwise  Jj -direction  giving  the  following  relation: 


-  1  / ( J  Re)  (  Cjiji  +  ^ )  =  -1/(J  Re)  !  +  S  /J 


(3.15) 


where 


l5=(Cx35a  +  5y35B  +  5z35c 


( 
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Lfl=(nx^r)^  +  Dy3qB  +  'r1z3qC)/J 
L^  =  ( +  i^yS^B  ■*■  5z3tC)/J 

Lnfi  =  nx3n  {^+vt)nx3ri^J  +r1y^r|{^+vt)r1y3r|D}  +T1z^ri  {^+vt)T1z9n^^ 

^CC  =  ^x3?  {( +Cy  3^  {( 1+',t)^y^c^>^  +^z^c  t  ^+vt)^z^C^l 

^ti5  =  t1x3ti {O+^t) Cx3 +rly3^  {( 1+^t^y^^l  +T1z3  p  {( 1+\>t)  CzS^D} 

^x3^{(^vt)rix3nD}  +  £y 3 £  { ( 1+vt) Dy 3 pD}  +^z3^  {( 1+Vt)hz3pD} 

where  J  is  the  Jacobian  of  the  transformation  given  by: 

J  =  £x(  hyCz-rlz£y )  ~  Dx( £y£z-^z£y ^  +  £x(£yDz“£zDy ) 

(3.16) 

and  vt  is  the  eddy  viscosity  described  in  chapter  IV.  Also  notice  that 
the  cross  derivative  term  is  held  explicitly  since  the  Linearized 
Block  Implicit  (LBI)  solution  procedure  of  Briley  and  McDonald  (1980) 
cannot  treat  such  terms  implicitly. 

The  vectors  in  equation  3.15  are: 
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:Cz(Pi+in*1-Pi-in)  +  n2 z-2ftv 

Equation  3.15  is  discretized  using  a  two  point  backward  difference 
scheme  for  streamwise  derivatives  and  a  three  point  central  difference 
scheme  for  all  transverse  derivatives.  A  linearized  form  of  the 
governing  equation  is  obtained  by  using  truncated  Taylor  series  about 
i-1  as  follows: 


A £  =  A^_ i  +  9qA  Aq 

(3.17) 

where  3qA  is  the  Jacobian  denoted  by  Aj  and  Aq  is  the  streamwise  change 
in  the  dependent  vector. 

T 

Aq  =  [ Ap,Au, Av, Aw] 

(3.18) 


Aq  =  qi  -  qi-l 


(3.19) 


When  equation  3.17  is  applied  to  all  the  vectors  in  equation  3.15  except 
S,  the  standard  delta  form  of  the  discretized  equation  results  giving: 

(  •  )  Aqi  =  R 

(3.20) 


where 


=  convection  terms  in  ^-direction 
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=  convection  and  diffusion  terms  in  q-direction 
=  convection  and  diffusion  terms  in  ^-direction 
R  =  terms  from  Euler  implicit  differencing  and  cross 
derivative  diffusion  terms 

The  operator  on  the  left  side  of  equation  3.20  is  split  into  two 
components  which  yields  block  tridiagonal  matrices  which  can  be  solved 
easily.  Solving  for  Aq  in  this  way  is  the  basic  LBI  scheme  of  Briley 
and  McDonald  (1980)  which  is  essentially  ADI  for  a  vector  system.  The 
method  is  second  order  accurate  in  the  transverse  plane  but  first  order 
accurate  in  the  streamwise  direction.  Full  second  order  accuracy  can  be 
obtained  using  three-point  backward  differences  for  the  streamwise 
gradients.  Since  truncated  Taylor  series  are  used  to  discretize  the 
governing  equation  and  no  averaging  is  used,  the  discretization  is 
consistent  by  definition.  The  method  does,  however,  introduce  an 
inconsistency  in  the  form  of  a  splitting  error. 

3.6  New  Solution  Procedure 


The  discretized  Navier-Stokes  equation  appears  to  be  second  order 
accurate.  However,  when  using  the  LBI  scheme  of  Briley  and  McDonald 
(1980),  a  splitting  error  arises  which  may  reduce  the  overall  accuracy 
of  the  method.  When  the  LBI  scheme  is  applied  to  a  space-marching 
method,  the  splitting  error  not  only  reduces  accuracy  but  may  destroy 
the  solution  due  to  an  inconsistency  in  the  splitting  error.  This  error 
can  be  examined  by  using  the  PL  method.  Equation  2.2  can  be  written  in 
frozen  coefficient  form  as: 


( 
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(Aj-  -  ACBj3n‘  +  A£Cj3c‘)Aq  =  r 


(3.21) 


The  diffusion  terms  are  neglected  to  simplify  the  analysis.  Equation 
3.21  is  generated  using  an  Euler  implicit  differencing  technique  and  the 
operator  on  the  left-hand  side  must  be  split  to  obtain  matrices  which 
can  be  inverted  easily.  The  split  can  be  carried  out  in  one  of  the 
following  two  ways. 


System  one 


first  solve  (Aj*  +  A£Bj3n»)Aq  =R 


then  solve  ( A j •  +  A^Cj3^  •  )Aq=AjAq 


(3.22) 


which  is  the  system  used  by  Pouagare  and  Lakshminarayana  (1986),  or 


System  two 


first  solve  ( A  j  •  +  A^Cj3^»)Aq  =R 


then  solve  (Aj*  +  A^BjS^ • )Aq=AjAq 


(3.23) 


When  system  one  and  two  are  recombined,  their  respective  splitting 


errors  are: 


System  one 


The  splitting  error  =  A£2  (Bj  Aj-1Cj)  3^*  3^»  Aq 


(3.24) 


System  two 


The  splitting  error  =  A?2  (Cj  Aj-1Bj)  3^*  3^*  Aq 


t 


77 


(3.25) 

The  splitting  errors  in  equations  3.24  and  3.25  are  equivalent  if  and 
only  if  the  matrices  Bj  and  Cj  are  commutative  matrices  which  they  are 
not  in  general.  This  was  recognized  by  Briley  and  McDonald  (1980)  but 
was  not  considered  to  be  a  problem  since,  for  time-marching  methods,  the 
splitting  errors  are  multiplied  by  a  time  change  in  the  dependent  vector 
which  tends  to  zero  for  steady  flows.  For  parabolic-marching  methods, 
the  splitting  error  is  multiplied  by  the  streamwise  change  in  the 
dependent  vector  which  seldom  goes  to  zero.  Thus,  the  splitting  error 
goes  to  zero  only  if  the  streamwise  grid  is  refined.  Refining  the  grid 
by  one  order  of  magnitude  will  give  a  two  order  reduction  in  the 
splitting  error.  Recall  though,  that  adding  ten  times  more  streamwise 
grid  lines  could  increase  the  number  of  global  iterations  of  a 
parabolic-marching  scheme  by  as  much  as  thirty  times.  Also,  a  more 
refined  grid,  while  giving  increased  accuracy  will  require  much  more 
computer  time  per  global  iteration.  The  magnitude  of  this  error  is 
unknown  and  so  should  be  investigated. 

The  matrix  multiplication  found  in  equations  3.24  and  3.25  can  be 


performed  for  the 

PL 

method. 

For  system  one 

"  0 

0 

w/u 

0 

Bj  Aj-’Cj  = 

0 

0 

w 

V 

0 

w/o 

2wv/u 

-u/o 

v/u 

0 

w2/u 

vw/u 

(3.26) 
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For  system  two 

"  0 

0 

0 

v/u 

Cj  Aj^Bj 

0 

0 

w 

V 

= 

w/u 

0 

wv/u 

v2/u 

.  0 

v/o 

-u/o 

2vw/u 

(3.27) 

Recall  that  the 

order  of 

the  equations  is  continuity,  x-,  y-,  z- 

momentum;  thus, 

the 

error 

is  the 

same  for 

the  x-momentum  equation  and 

very  similar  for  the  y-  and  2-momentum  equations.  However,  the  error  in 
the  continuity  equation  is  much  different  between  the  two  systems. 

With  the  elements  of  Aq  being  roughly  the  same,  if  v  is  much  larger 
than  w  as  would  be  the  case  for  a  duct  curving  in  the  x-y  plane,  then 
one  expects  system  one  to  give  a  smaller  error  in  the  continuity 
equation  than  system  two  and  thus  better  results.  This  observation  will 
be  born  out  in  section  3-6.3.  One  other  problem  that  should  be  noted  is 
that,  since  a  forward-marching  procedure  is  under  consideration,  the 
error  will  build  up  with  each  marching  step.  There  is  no  mechanism  to 
reduce  the  error  as  the  computation  progresses. 

3.6.2  New  Solution  Method 

The  error  can  be  significantly  reduced  by  iterating  equation  3.21  at 
each  streamwise  station.  For  convenience  the  equation  under 
consideration,  i.e.,  equation  3.20,  is  repeated  here. 

(  L^*  +Ln*  +L.£ •  )  Aqi  =  R 

(3.28) 


An  iterative  form  of  3.28  can  be  written  as: 


79 


(  L^*  +Lrj*  +L^*  )  Asim  =  R 
R'=  R  -  (  Lj-  •  +Ln*  +L^  *  )  Aqim_1 

(3.29) 


where  Asim=Aqjm-Aqjn1-1  with  m  as  an  iteration  index  and  R*  being  the 
solution  to  the  Navier-Stokes  equation  within  the  formal  accuracy  of  the 
discretization . 

The  operator  on  the  left  side  of  equation  3.29  can  now  be  split  as 
in  the  LBI  scheme,  however,  the  splitting  error  is  reduced  with  each 
iteration  at  streamwise  station  i  since  As^™  tends  to  2ero.  To  improve 
convergence,  Aqim  is  updated  using  the  following  relationship. 


Aqim  =  Aqim-^  +  w  Asim 

(3.30) 

where  w  is  an  underrelaxation  parameter  which  usually  is  set  to  a  value 
of  0.7. 

The  solution  procedure  for  equation  3.29  is  as  follows.  At 
streamwise  station  i,  the  operator  in  equation  3-29  is  split  using 
either  system  one  or  two  using  R  as  determined  from  the  previous 
streamwise  station  and  the  value  of  Aqim_1;  for  the  first  iteration 
this  is  the  converged  value  of  Aq  at  i-1.  With  the  right-hand  side 
known,  equation  3.29  is  solved  for  Asim  which  is  used  to  update  Aqj01 
using  equation  3.30.  The  procedure  is  repeated  until  As  drops  at  least 
two  orders  of  magnitude.  Such  convergence  is  usually  achieved  in  five 
iterations.  After  convergence,  the  procedure  is  advanced  to  the  next 
streamwise  station  and  repeated  until  the  entire  domain  is  computed.  In 
short,  the  splitting  error  is  iterated  out  of  the  solution  at  each 
streamwise  station  and  thus  the  discretized  system  retains  its  formal 


« 


4 
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second  order  accuracy. 

3.6.3  Error  in  Curved  Duct  Flow 

In  order  to  examine  the  effect  of  the  splitting  error  and  the 
ability  of  the  new  solution  method  to  reduce  it,  the  laminar  flow  in 
Taylor  et  al.'s  (1981)  strongly  curved  duct  was  computed.  The  geometry 
is  given  in  Figure  3.7.  The  radius  of  curvature  Rc  is  2. 3D  and  the  flow 
Reynolds  number  based  on  D  and  the  bulk  inlet  velocity  is  790.  The 
inlet  velocity  field  was  interpolated  from  the  available  experimental 
data  and  a  radial  pressure  gradient  based  on  simple  radial  equilibrium 
was  used  as  the  assumed  pressure  for  the  PL  method.  The  marching 
procedure  was  started  at  the  inlet  plane  <J>  =  0°  and  progressed  to  4>  = 
78°.  The  computational  grid  included  21  points  in  the  radial 
direction,  11  points  from  the  lower  wall  to  the  symmetry  plane  and  A<t>  = 
2°.  Recall  that  the  major  motivation  for  the  new  solution  method  is  to 
maintain  a  high  degree  of  accuracy  with  a  coarse  grid  for  improved 
convergence  of  the  parabolic-marching  method.  The  afore  mentioned 
splitting  error  is  small  for  fine  grids.  Indeed,  the  prediction  of  this 
flow  by  Pouagare  and  Lakshminayana  (1986),  reprinted  in  Figure  3.8,  is 
very  accurate  for  a  grid  of  41  by  21  in  the  cross  plane  and  A$=  0.5°, 
nearly  15  times  as  many  points  as  the  present  computation.  The  method 
used  by  Pouagare  and  Lakshminarayana  is  equation  3.28  using  the  system 
one  splitting  procedure.  This  method  was  used  as  well  as  that  using  the 
splitting  procedure  of  system  two  on  the  coarse  grid  and  compared  to  the 
results  obtained  by  the  new  iterative  solution  method. 

Computed  streamwise  and  transverse  velocity  profiles  at  $  =  77.5° 
are  presented  in  Figures  3.9  and  3.10,  respectively.  Careful 


Figure  3.8  Fine  Grid  Results  of  Pouagare  and  Lakshminarayana  (1986) 
for  the  Strongly  Curved  Duct  Flow 


Solution  method - -  System 

Taylor  (1981)  -  System 


Computed  Streamwise  Velocity  Profiles  for  the  Strongly 
Curved  Duct  Flow,  PL  method 


Curved  Duct  Flow,  PL  method 
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examination  of  the  results  shows  the  pronounced  effect  of  the  splitting 
error.  The  results  using  system  one  agree  reasonably  well  with  the 
experimental  data  as  does  those  of  the  new  solution  method.  The  results 
from  using  system  two  are  very  poor  as  was  predicted  from  the  analysis 
in  the  previous  section.  The  difference  between  the  two  systems  is 
quite  dramatic  and  is  due  soley  to  the  effect  of  the  splitting  error. 
The  new  solution  method  produces  results  closer  to  the  experimental  data 
especially  in  the  secondary  flow  predictions.  Note  the  large  deviation 
from  the  data  at  y/D  =  0.9  in  Figure  3.10  for  system  one  while  the 
results  of  the  new  solution  procedure  compare  very  well. 

The  residual  in  the  continuity  equation  is  plotted  in  Figure  3.11 
verses  the  streamwise  station  for  the  three  methods  of  splitting  the 
operator.  Notice  that  the  residual  for  system  one  is  much  less  than  for 
system  two  as  anticipated,  however,  the  error  is  clearly  unbounded  as 
the  computation  is  carried  farther  downstream.  On  the  other  hand,  the 
error  in  the  continuity  equation  using  the  new  iterative  solution  method 
is  much  lower  and  remains  low.  It  is  evident  from  this  result  that  the 
extra  computational  effort  required  by  the  new  solution  method  is 
Justified. 

The  convergence  history  of  the  iterative  LBI  scheme  at  a  typical 
streamwise  station  for  this  flow  is  given  in  Figure  3.12.  Strong 
convergence  is  exhibited  for  all  variables.  It  should  be  noted  that 
system  one  and  two  were  used  in  alternative  iterations  and  is  the  reason 
for  the  small  spikes  in  the  convergence  history.  Clearly,  one  splitting 
system  is  better  suited  to  this  flow  but  the  iterative  scheme  still 
converges . 


Figure  3.11  Computed  Residual  in  the  Continuity  Equation 
for  the  Three  Solution  Methods 
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3.7  Boundary  Conditions 

For  the  internal  turbomachinery  problems  for  which  the  NPM  method 
was  developed,  only  four  different  boundary  conditions  are  considered. 
First,  the  no-slip  boundary  condition  is  the  proper  boundary  condition 
for  viscous  flow  near  a  solid  boundary.  For  such  a  condition,  all 
velocities  are  set  to  zero  on  the  boundary  with  the  normal  gradient  of 
the  pressure  set  to  zero  following  boundary  layer  theory.  A  two  point 
backward  difference  scheme  is  used  for  the  pressure  gradient  condition 
in  order  to  couple  the  odd  and  even  points  more  strongly  (see  Pouagare 
and  Lakshminarayana  (1986)).  This  condition  is  always  used  for  laminar 
flows.  Second,  for  turbulent  flows,  the  no-slip  condition  may  also  be 
used.  However,  under  special  circumstances,  namely  high  turhulance 
intensities  with  a  coarse  grid,  a  turbulent  slip  velocity  must  be  used. 
This  slip  velocity  is  based  on  the  log  law-of-the-wall  and  is  determined 
as  follows: 


Gwall  =  Qp  ~  “»/< 


(3.3D 


where  Q  is  the  total  velocity,  ic  is  the  von  Karman  constant,  and  u#  is 
the  friction  velocity  determined  from  the  following  relationship: 


Qp/u«  =  5.8  log  y+  +  5.0 


(3.32) 


where  y+  =  Re  u*  np 

(3.33) 

np  is  the  normal  distance  from  the  boundary  to  the  first  grid  point. 
This  condition  is  carefully  applied  so  that  the  resultant  velocity 
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component  is  tangent  to  the  boundary  surface. 

Third,  on  a  symmetry  plane,  the  normal  gradient  of  the  velocities 
tangent  to  the  plane,  one  of  which  is  the  streamwise  component,  is  set 
to  zero.  The  velocity  component  normal  to  the  plane  of  symmetry  is  also 
set  to  zero.  The  normal  pressure  gradient  is  set  to  zero  across  the 
symmetry  boundary.  A  two  point  backward  difference  scheme  is  employed 
for  these  gradient  conditions  both  implicitly  and  explicitly. 

Finally,  for  most  of  the  turbomachinery  flows  investigated  here, 
there  are  substantial  regions  of  the  flow  which  are  spacially  periodic. 
Thus,  a  periodic  boundary  condition  is  used  upstream  and  downstream  of  a 
blade  row.  This  condition  can  be  enforced  in  several  ways.  First,  the 
condition  can  be  enforced  implicitly.  This  requires  the  solution  of  a 
periodic  matrix  which  can  be  solved  using  a  recursive  formula.  The 
matrix  has  the  form: 

be  a 

a  b  c 

c  a  b_ 

Without  artificial  dissipation,  central  difference  schemes  notoriously 
produce  uncoupled  solutions.  The  solution  of  the  equations  with 
implicit  periodic  boundary  conditions  exacerbates  the  odd  even 
uncoupling. 

Second,  the  periodic  condition  can  be  applied  explicitly.  The 
computational  domain  can  be  extended  to  include  the  periodic  line  on 
both  sides  of  the  domain.  The  value  of  q  computed  at  the  previous 
iteration  in  the  transverse  plane  is  used  as  boundary  conditions  outside 


g 
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the  periodic  lines.  Thus,  two  different  values  of  q  might  be  computed 
on  the  periodic  line,  an  upper  and  a  lower  value.  To  maintain 
periodicity,  the  upper  and  lower  values  are  averaged  to  give  the  value 
of  q  on  the  periodic  line.  This  slows  convergence  of  the  iterative 
scheme  in  the  transverse  plane  but  periodicity  must  be  enforced. 
Uncoupling  of  the  odd  and  even  points  is  not  as  severe  when  the  explicit 
condition  is  applied.  It  should  be  noted  that  for  coupling  purposes, 
the  normal  pressure  gradient  equals  zero  condition  is  always  applied 
implicitly  regardless  of  which  boundary  condition  is  used.  Experience 
has  shown  that  there  is  no  great  loss  in  accuracy,  even  in  periodic 
regions,  when  this  condition  is  used,  however,  the  coupling  of  the  odd 
and  even  points  is  greatly  enhanced. 

At  the  inlet  plane,  only  the  velocity  components  are  needed  as  inlet 
boundary  conditions.  This  inlet  field  can  be  generated  from  available 
experimental  data.  Downstream,  only  the  pressure  is  used  as  a  boundary 
condition  and  so  needs  to  be  specified.  Since  the  flow  is 
incompressible,  the  level  of  pressure  is  not  important  but  the 
transverse  pressure  gradients  downstream  should  be  correct,  i.e.,  match 
the  experimental  data  or  analysis. 


The  NPM  algorithm  must  be  calibrated  for  convergence  and  accuracy  on 
a  simple  laminar  flow  before  one  can  have  faith  in  its  ability  to 
predict  the  complex,  turbulent  flow  in  a  turbomachine.  The  most  simple 
yet  non-trivial  laminar  flow  that  can  be  computed  is  the  developing  flow 
in  a  straight  duct  with  a  square  cross-section.  Once  a  successful 
computation  of  this  flow  has  been  achieved,  the  laminar  flow  in  an  S- 
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shaped  duct  of  square  cross-section  is  computed  since  this  flow  is 
characterized  by  large  streamwise  variations  in  the  transverse  pressure 
gradients  as  well  as  large  secondary  flows.  These  characteristics  are 
also  found  in  most  turbomachinery  flows. 

3.8.1  Developing  Laminar  Flow  in  a  Straight  Duct 

For  this  case  the  Reynolds  number  based  on  the  duct  width  was  set  to 
50.  The  domain  was  extended  to  6  duct  widths  downstream  where  fully 
developed  flow  should  be  present.  A  grid  density  of  21x21  points  was 
used  in  the  cross-plane  while  25  points  were  used  in  the  streamwise 
direction.  Exponential  stretching  was  used  to  pack  grid  points  near  the 
entrance  region.  A  uniform  inlet  profile  was  used  to  start  the  marching 
process  and  a  unifrom  pressure  field  was  used  as  the  downstream  pressure 
boundary  condition.  For  these  conditions,  the  optimum  value  for  a  is 
computed  by  equation  3.14b  is  a=0.2.  The  value  of  a  was  varied  from  0 
to  0.4  in  order  to  study  its  effect  on  the  rate  of  convergence.  Figure 
3.13  shows  the  convergence  history  for  a  =  0,  0.01,  0.2  and  0.4.  It  is 
evident  that  for  the  initial  passes,  convergence  is  faster  for  a  not 
equal  to  zero  as  the  analysis  professed.  Overall,  the  fastest 
convergence  was  achieved  with  a=0.0.  This  is  not  suprising  since  the 
normal  pressure  gradients  are  essentially  zero  and  the  satisfaction  of 
the  global  mass  flow  constraint,  through  continuity,  is  the  driving 
mechanism  in  this  particular  flow.  Therefore,  with  a=0.0,  convergence 
should  be  achieved  in  as  many  global  iterations  as  there  are  streamwise 
points  which  is  indeed  the  case.  Still,  with  a  not  equal  to  zero, 
convergence  is  achieved  (initially  at  a  faster  rate)  but  clearly  a=0.2 
is  not  the  optimum  value.  This  gives  credence  to  the  observation  that, 


while  an  optimum  value  for  a  apparently  can  be  computed  for  one  set  of 
wave  numbers  in  the  spectrum,  it  may  not  be  the  optimum  value  for  all 
wave  numbers. 

The  computed  centerline  coefficient  of  pressure  and  the  centerline 
velocity  are  compared  to  the  analytical  results  (see  White  (1974))  in 
Figures  3.14  and  3.15,  respectively.  One  can  see  that  the  computed 
coefficient  of  pressure  compares  favorably  with  the  analytical  results 
except  very  near  the  duct  entrance.  The  computed  centerline  velocity 
compares  well  with  the  analytical  results  only  far  from  the  entrance. 

In  the  entrance  region,  the  velocity  profiles  tend  to  bulge  near  the 
walls  of  the  duct  reducing  the  momentum  in  the  center  of  the  duct.  This 
is  most  likely  due  to  the  development  of  spurious  transverse  pressure 
gradients  which  suppress  the  diffusion  of  momentum  from  near  the  wall  to 
the  center  of  the  duct.  This  did  not  occur  when  the  PL  technique  was 
used  by  Pouagare  and  Lakshminarayana  (1986).  This  is  because  in  the  PL 
technique,  the  small  value  of  o  essentially  removes  the  effect  of  the 
transverse  pressure  gradients  from  the  computation.  These  spurious 
transverse  pressure  gradients  are  not  large  and  diminish  to  zero  past 
the  entrance  region.  Grid  refinement  in  the  entrance  region  should  be 
all  that  is  necessary  to  reduce  this  problem.  In  the  fully  developed 
flow  region,  the  exact  centerline  velocity  is  computed  and  the  computed 
velocity  profiles  compare  very  well  to  the  analytical  results  in  Figure 
3.16.  These  results  did  not  change  for  the  various  values  of  a.  This 
proves  that  at  convergence,  the  continuity  equation  is  satisfied  within 
the  accuracy  of  the  discretization  everywhere  and  the  global  mass  flow 


is  conserved. 


x/Re/D 


Figure  3.15  Computed  Centerline  Velocity  for  the  Straight  Duct 
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3.8.2  Laminar  Flow  in  an  S-Shaped  Duct 

The  entry  flow  in  a  straight  duct,  while  not  trivial,  is  not 
characterized  by  a  streamwise  variation  in  the  transverse  pressure 
gradients  as  is  a  typical  turbomachinery  flow.  Thus,  the  laminar  flow 
in  an  S-shaped  duct  was  computed  to  determine  the  NPM  method's  ability 
to  resolve  the  large  variations  in  the  transverse  presure  gradients  as 
well  as  the  strong  secondary  flows  present.  The  flow  was  measured  by 
Taylor  et  al.  (1982)  using  a  laser  Doppler  velocimeter.  The  geometry  is 
given  in  Figure  2.9  along  with  the  measurement  locations.  The  radius  of 
curvature  was  7D  to  the  center  of  the  duct.  The  Reynolds  number  based 
on  the  bulk  inlet  velocity  Ug  and  the  duct  width  D  was  790.  The 
computational  domain  extended  from  Xh  =  -2.5  to  8.0  and  from  zero  to  D 
in  both  the  Y  and  Z  directions.  The  comptuational  grid  consisted  of  47 
points  in  the  streamwise  direction  spaced  evenly  and  a  25  by  21  grid  was 
used  in  the  cross-plane  with  exponential  stretching  included  to  pack  the 
grid  points  near  the  walls.  The  inlet  conditions  were  determined  by 
computing  the  developing  flow  in  the  entrance  region  up  to  Xh=-2.5  using 
the  PL  technique. 

Recall  that  for  the  developing  flow  in  a  straight  duct,  the  best 
convergence  was  achieved  with  a=0.0.  With  a=0.0,  the  governing  equation 
is  essentially  unmodified  and  the  method  is  no  different  than  the  PL 
method  with  o=-1.0.  However,  with  no  assumed  pressure  distribution, 
the  computation  could  not  progress  to  the  end  of  the  domain  successfully 
with  a=0.0.  With  a  set  to  5.0,  however,  the  computation  was  successful. 
Figure  3.17  shows  the  convergence  history  for  a=5.0.  Convergence  to  two 
orders  of  magnitude  were  achieved  in  140  global  iterations.  Total  CPU 
time  on  the  IBM  3090-200  at  The  Pennsylvania  State  University  was  1.87 
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hours  or  0.00195  seconds  per  point  per  global  iteration. 

Figure  3.18  shows  the  computed  pressure  distribution  for  the  S-duct. 
No  experimental  data  is  available  for  the  Cp  distribution  for  the 
laminar  flow  but  the  computed  head  loss  is  very  close  to  that  for  the 
developing  flow  in  a  straight  duct  with  the  same  inlet  conditions.  The 
results  are  certainly  qualitatively  correct  with  a  negative  transverse 
pressure  gradient  in  the  upward  curving  region  switching  to  a  positive 
gradient  as  the  duct  turns  back.  Indeed,  when  comparing  the  computed  Cp 
to  that  measured  for  turbulent  flow  (see  Figure  2.14),  the  fine  details 
of  the  distribution  agree  very  well  qualitatively. 

The  computed  velocity  profiles  are  compared  to  the  experimental  data 
in  Figures  3.19  through  3.26.  The  agreement  between  the  computation  and 
the  experimental  data  is  excellent  in  all  regions  of  the  flow.  At 
station  1,  Just  upstream  of  the  initial  bend,  the  computed  velocity 
field  is  no  longer  symmetric  and  a  small  secondary  velocity  exists  owing 
to  the  transverse  pressure  gradient  effect.  Since  streamwise 
ellipticity  was  reduced  to  first  order  in  the  NPM  method,  it  is  apparent 
that  the  pressure  ellipticity  is  sufficient  for  an  accurate  prediction 
for  this  type  of  flow.  At  station  2,  a  strong  secondary  flow  is 
developing  due  to  the  influence  of  the  centrifugal  force  and  the  end 
wall  boundary  layer.  This  is  typical  in  turbomachinery  blade  rows. 
Again,  agreement  with  the  experimental  data  is  excellent  at  all 
measurement  locations.  At  station  4,  the  reversal  of  the  secondary  flow 
due  to  the  change  in  sign  of  the  transverse  pressure  gradient  has  begun. 
Also  apparent  is  the  bulge  in  the  streamwise  velocity  profile  at  Y/D=0.7 
indicating  the  center  of  the  secondary  vortex.  Again  the  method 
accurately  predicts  all  these  flow  phenomena.  Finally  at  station  5,  the 
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Figure  3.18  Computed  Coefficient  of  Pressure  for  the  S-Duct 
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secondary  flow  reversal  is  complete  and  the  influence  of  the  straight 
portion  of  the  duct  outlet,  through  the  streamwise  pressure  gradient, 
has  begun  to  straighten  out  the  flow  at  Y/D=0.9.  At  this  location,  the 
streamwise  profile  is  not  well  predicted  near  the  symmetry  plane.  It 
should  be  noted  that  the  grid  is  sparse  in  this  region  and  may  have  led 
to  the  poor  prediction.  This  exemplifies  one  of  the  major  obstacles 
hindering  a  true  prediction  of  the  flow.  If  experimental  data  were  not 
present,  the  secondary  vortex  that  generates  the  odd  profile  at  Y/D=0.9 
would  not  be  anticipated  and  a  fine  grid  would  normally  not  be  placed  in 
that  region. 

Overall,  the  agreement  with  the  experimental  data  for  the  S-duct  was 
excellent.  This  case  proved  that  multiple  passes  of  the  partially 
parabolized  Navier-Stokes  equation  as  a  solution  procedure  is  capable  of 
transmitting  downstream  pressure  effects  upstream  albeit  slowly  and  that 
the  NPM  method  converges  even  when  all  the  pressure  gradients  are 
substantial.  Also,  the  case  proved  that  the  NPM  method  is  superior  to 
multiple  pass  methods  which  only  forward  difference  the  streamwise 
pressure  gradient  with  no  other  modifications  since  the  flow  could  not 
be  computed  with  a=0.0. 


3.9  Closing  Remarks 

In  the  most  general  terms  the  NPM  method  is  a  line  relaxation  of 
Chorin's  method  with  no  unsteady  velocity  terms.  With  the  inclusion  of 
the  pressure  term  in  the  continuity  equation,  it  has  been  shown  that 
convergence  is  improved  over  a  partially  parabolic  system  with  no 
modifications.  Also,  the  extra  term  tends  to  improve  the  coupling  of 
the  odd  and  even  points.  In  addition,  with  a  non-zero  element  in  the 
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1,1  position  of  Aj,  diagonal  dominance  is  assured  and  the  amount  of 
pivoting  required  in  the  solution  of  the  block  tridiagonal  matrices  is 
reduced  considerably  with  a  saving  of  computation  time.  It  should  be 
noted  that  no  implicit  or  explicit  smoothing  is  required  for  good 
results.  Due  to  the  one-sided  differencing  of  the  pressure  gradient 
boundary  condition  and  the  one-sided  differencing  of  all  streamwise 
derivatives,  there  is  a  small  amount  of  artificial  dissipation 
introduced. 

The  method  has  been  shown  to  accurately  introduce  the  streamwise 
ellipticity  due  to  the  pressure  as  seen  in  the  S-duct  results.  This  is 
the  most  important  feature  of  the  method  since  strong  pressure 
ellipticity  is  the  reason  one  cannot  use  the  PL  or  MPL  method  with 
success.  The  method  has  been  used  to  compute  laminar  flows  with 
excellent  results  and  so  it  should  be  able  to  compute  complex  turbulent 
flows  with  success  with  the  inclusion  of  turbulence  models. 

Due  to  the  parabolizing  assumptions  made  in  the  governing  equation, 
only  H-type  computational  grids  can  be  used.  A  description  of  this  type 
of  grid  and  an  algorithm  for  generating  them  is  described  in  the 
appendix.  Also,  the  flow  near  blunt  leading  edges  cannot  be  computed 
accurately  due  to  the  lack  of  streamwise  diffusion  in  the  governing 
equation. 


Ill 


CHAPTER  4 

TURBULENCE  CLOSURE  SCHEMES 

4.1  Introduction 

Direct  simulation  of  the  turbulence  by  solving  the  unsteady  Navier- 
Stokes  equation  at  the  Kolmogorov  scales  is  a  monumental  task  for 
complex  turbomachinery  flows.  Thus,  a  simplified  form  of  the  equation, 
namely  the  Reynolds  averaged  Navier-Stokes  equation,  is  used  which  is 
similar  to  equation  2.1  with  the  variables  being  time  averages.  The 
Reynolds  averaging  also  introduces  time  averaged  products  of  velocity 
fluctuations  called  Reynolds  stresses.  These  Reynolds  stresses  can  be 
interpreted  as  apparent  additional  shear  stresses  and  related,  using  the 
Boussinesq  approximation,  to  a  mean  velocity  gradient  with  some  constant 
of  proportionality.  This  constant  is  termed  the  eddy  viscosity  since  it 
mimics  the  kinematic  viscosity  but  is  associated  with  the  turbulence 
eddies.  With  this  model,  the  governing  equation  is  in  a  closed  form 
provided  the  eddy  viscosity  is  known.  Here  lies  the  difficulty  in 
modelling  turbulence.  The  eddy  viscosity  must  somehow  be  modelled  using 
our  knowledge  and  intuition  about  the  nature  of  turbulence. 

Various  models  have  been  proposed  with  varying  degrees  of 
simplification  and  empericisra.  The  most  simple  model  is  the  algebraic 
eddy  viscosity  model  which  relates  the  eddy  viscosity  to  some  length 
scale  and  a  mean  velocity  gradient.  The  two-equation  k-e  model  relates 
the  eddy  viscosity  to  the  turbulence  kinetic  energy  and  the  turbulence 
energy  dissipation.  Both  of  these  models  do  not  include  the  effects  of 


the  anisotropy  of  turbulence.  Unfortunately,  turbulence  in  general  is 
not  isotropic.  Flow  curvature,  rotation  and  bouyancy  all  introduce 
anisotopy  which  should  be  modelled  appropriately.  The  best  way  to  do 
this  is  to  use  the  full  Reynolds  stress  model  where  partial  differential 
equations  are  solved  for  each  Reynolds  stress.  Some  of  the  terms  in 
these  equations  must  be  modelled  themselves  which  tends  to  complicate 
matters  greatly.  A  compromise  between  the  k-e  model  and  the  full 
Reynolds  stress  model  is  the  algebraic  Reynolds  stress  model  or  ARSM. 
The  ARSM  uses  algebraic  equations  for  the  various  Reynolds  stresses  and 
includes  some  degree  of  anisotropy.  The  ARSM  along  with  the  algedraic 
eddy  viscosity  model  and  the  k-e  model  will  be  discussed. 


4.2  Algebraic  Eddy  Viscosity  Model 


From  a  purely  dimensional  argument,  the  eddy  viscosity  can  be 
related  to  some  turbulence  length  scale  4  and  a  mean  velocity  gradient. 
Thus,  the  equation  for  the  eddy  viscosity  is 


vt  =  42  3nQ 


(4.1) 


A  two  layer  algebraic  model  developed  by  Baldwin  and  Lomajc  (1978)  has 
been  used  here  for  all  shear  layers  where: 


»t  = 


vt  inner  0  <  n  <  nc 


^t  outer  nc  <  n  <  6 


(4.2). 


where  n  is  the  distance  normal  to  the  boundary,  5  is  the  edge  of  the 
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shear  layer  and  nc  is  the  point  where  the  eddy  viscosity  is  switched 
from  the  inner  to  the  outer  formula.  The  length  scale  for  the  inner 
formulation  includes  the  van  Driest  damping  function  and  the  total 
vorticity  replaces  the  mean  velocity  gradient  in  equation  4.1.  Thus, 
the  non-dimensional  inner  eddy  viscosity  is  written  as: 


vt  inner  =  C<  n  ( 1-exp( -y"7A+) ) ]2 | w  |  Re 


(4.3) 


where 

b> 


/ 


(3yU-3xv)^  +  {3zv-3yw)^  +  (3xw-3zu)‘ 


(4.4) 


and 

y+  =  u*  n  Re 
k  =  0.4 
A+  =  26 

u*  =  friction  velocity  defined  in  equation  3.32 

The  outer  formulation  is  based  on  several  emperical  constants,  the 
Klebenov  intermittancy  factor,  and  the  total  mean  vorticity. 


vt  outer  =  K  CCp  Fwake  ^kleb 


(4.5) 


where 


Fwaice  =  minimum  of 


nmax  Fmax 


Cwk  ^max  Qd  /Fmax  » 
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nmax  is  the  distance  where  F  is  a  maximum  from: 
F  =  <  n  (1  -  exp(-y+/A+)) |u  j 
^kleb  =  [  1  +  5.5(  C^ieb  n/nmax)^]  ^ 


and 

Qd  =  Qmax  ~  Qmin 

The  following  emperical  constants  are  applied. 

CCp  =  1.6 

Ckleb  =  0-3 

Cwk  =  0.25 

K  =  0.0168 

In  determining  the  crossover  location  nc,  it  has  been  found  to  be 
more  computationally  efficient  to  compute  both  the  inner  and  outer 
components  of  the  eddy  viscosity  and  using  the  smaller  of  the  two.  In 
three  dimensions,  the  above  formulations  were  used  for  the  nearest 
boundary.  It  should  also  be  noted  that  the  distribution  of  F  has  two 
peaks.  Experience  has  shown  that  taking  Fmax  to  be  the  second  peak 
yields  the  best  results.  For  wakes,  only  the  outer  formulation  of 
is  used. 

The  advantage  of  the  Baldwin  and  Lomax  model  over  other  algebraic 
eddy  viscosity  models  is  that  the  edge  of  the  shear  layer  need  not  be 
computed  directly.  Rather,  the  shear  layer  edge  is  determined  from  the 
vorticity.  This  is  important  for  internal  flows  since  the  freestream  or 
the  inviscid  core  is  not  uniform  but  characterized  by  some  non-zero 
gradient.  Thus,  using  the  traditional  determination  of  the  shear  layer 


edge,  i.e.,  at  Q  =  .99  Qfreestream*  becomes  difficult. 

This  model  can  be  used  for  a  wide  variety  of  flows  with  very  good 
results.  Difficulties  arise  in  large  interaction  regions  ,  e.g., 
separated  flows,  and  highly  three-dimensional  flows.  The  algebraic 
model  contains  no  'history'  effects  of  the  turbulence  and  is  based  on 
two-dimensional  shear  layers  only.  In  order  to  include  more  physics  in 
the  model,  a  two-equation  k-e  turbulence  model  can  be  used. 


4 . 3  Two-equation  model 


The  greatest  difficulty  with  the  above  algebraic  model  is  that  the 
turbulence  length  scale  must  be  determined  in  an  ad  hoc  fashion.  It  is 
desirable  to  determine  the  turbulence  length  scale  and  velocity  scale 
from  partial  differential  equations  derived  from  physical 
considerations.  The  turbulance  velocity  scale  can  be  found  from  the 
scalar  turbulence  kinetic  energy.  An  equation  for  the  turbulence  length 
scale  does  not  exist  but  the  turbulence  dissipation  in  combination  with 
the  turbulence  kinetic  energy  does  yield  a  length  scale  and  a  partial 
differential  equation  can  be  written  for  the  turbulence  dissipation. 
From  dimensional  analysis,  the  eddy  viscosity  can  be  related  to  the 
turbulence  kinetic  energy  and  dissipation. 


vt  *  Vr/z 


in  non-dimensional  form,  equation  4.6  is  written  as: 


(4.6) 


Vr  =  Cn  k Ve  Re 
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The  equation  for  the  turbulence  kinetic  energy  can  be  derived  from  the 
Reynolds  stress  equation  and  has  the  following  form: 

9](Uk  +  dyVk  +  9xwk  —  1/Re  [  9x(vk3xk)  ^  ^y( ^k^yk) 

^k^z^  ]  =  Pj(/Re  -  e 

(4.8) 


where 

vk  =  1  +  vt/ok. 

The  turbulence  dissipation  equation  can  be  derived  by  taking  the 
derivative  of  the  unsteady  Navier-Stokes  equation  with  respect  to  the 
coordinates,  time  average,  then  subtract  the  mean  flow  equation  and 
multiply  the  result  by  the  unsteady  strain  rate  and  the  kinematic 
viscosity.  Taking  the  time  average  gives  an  equation  for  the  isotropic 
turbulence  dissipation.  The  equation  includes  several  terms  which  must 
be  modelled.  The  final  form  of  the  equation  has  the  following  form: 

3xue  +  3yV£  +  dgW£  •  1 /Re  [  9x(ve9xe)  ^  9y(vE9y£) 

+  Ce ie/k  Pk  ~  Ce2£  /k  Re  J  =  0 

(4.9) 


where 

ve  =  1  +  vt/oE. 

The  various  constants  were  set  to  match  experimental  data  for  free  and 
bounded  shear  layers  and  dissipation  of  turbulence  in  gradient  free 
flow.  The  standard  values  are  (see  Rodi  (1982)): 


CEl  =  1.44 


Ce2  =  1.92 
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0k  =  1 .0  0£  =  1 .3 

cu  =  0.09 

The  production  of  turbulence  is  given  by  the  following  simplified 
relationship: 

Pk  =  vttOn  Q>2  +  <3e  Q)2^ 

(4.10) 


The  k-e  equations  must  be  recast  into  generalized  coordinates  in  a 
manner  similar  to  the  mean  flow  equations.  Thus,  the  k-equation 
becomes : 

+  Lf|  +  —  1  / ( J  Re)  (  +  L^^)  =  —  ^/(J  Rs)  +  S*p/J 

(4.11) 


where 

L^=(^x3^uk  +  £y8rvk  +•  ^Z3^wk)/J 

Lq=( qx3  qUk  +  ny3nvk  +  hz9qWk)/<J 

L^  =  (  Cx3?uk  +  £y3^vk  +  k)/J 

4qq=hx^q  { vknx3nk^  +r1y3Ti  ivkHy3qk}  +rlz3q  { vkhz3r|l{} 

LCC=Cx3c  (vk^x3Ckl  +Sy3C  {vk^y3Ck}  +Cz3c  (vkCz3ck} 

l-TjC=nx3ii { vk5x3C^}  ^OySq {vjfCy 9^^}  -*-nz9q {vk^z3C^l 
+Cx3C  (vkT1x3hk}  +?y3C  (vkT1y3Tik}  +tz3C  {vkT'z3nkl 

St  =  Pk/Re  ~  e 


Similarly,  the  e-equation  becomes: 
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♦  Ln  +  LC  -  1/(J  Re)  (  Lnn  ♦  Lcc)  +  ST/J  =  -  1/(J  Re)  L nc 

(4.12) 


where 

L^=(^x3^ue  ♦  SyS^ve  +  5zHwe^J 
Lf|  =  (  hx  3  nu^-  +  Hy3r|V£  +  Dz3qw£)/J 
L^=(Cx^5u^  +  Cy3^ve  +  Cz3^we)/J 

^■riTi^^x^ri  {ver'x3nc  }  +r1y3ri  WeDyS^e }  +  r1z3 n  { v£T1z3nE  1 

LCC=£x3c {ve£x3ce}  +^y3C  {veCy3££  }  +  £z3C  NeCz9^} 

^tlC  =  rix^r|  (ve^x3^£ }  +Tly 9ri  f  ve5y3;;£ }  ■*■  Hz 3 rj  { ve£z^i;e } 

+i;x3^  { venx3  n£  J  +^y3^  {vcTiy3r|e}  +Cz^c  {ver>z^ne^ 

ST  =  ce1ePk/(k  Re)  -  Ce2£2/k 


The  k-e  equations  are  weakly  coupled  and  so  are  solved  in  an  uncoupled 
fashion  with  the  k  equation  solved  first.  Equations  4.11  and  4.12  are 
discretized  in  the  same  manner  as  the  mean  flow  equations  and  solved  in 
the  cross  plane  using  the  alternating  direction  implicit  (ADI) 
algorithm.  Ti  solution  of  the  equations  is  lagged  one  streamwise 
station  behind  the  mean  flow  equation  using  the  velocities  previously 
computed . 


4.3.1  Boundary  Conditions 

Since  the  k-£  equations  are  valid  only  in  regions  where  the 
turbulence  Reynolds  number  is  high,  wall  functions  are  used  to  give 
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boundary  conditions  for  the  turbulence  kinetic  energy  and  dissipation 
which  avoids  having  to  solve  the  k-e  equations  in  the  viscous  sublayer 
and  buffer  region  of  a  bounded  shear  layer.  The  boundary  conditions  are 
those  proposed  by  Rodi  (1982).  At  the  first  point  from  a  solid 
boundary,  the  following  values  of  k  and  e  are  applied  and  the  k-e 
equations  are  solved  in  the  rest  of  the  domain. 


(4.13) 


=  u*-5/ ( tc  nD) 


(4.14) 


where  u#  is  the  friction  velocity  defined  in  equation  3-32.  For  these 
wall  functions  to  apply,  the  first  point  from  the  solid  boundary  must  be 
no  closer  than  at  y+  =  20.  Although  the  use  of  these  wall  fuctions 
implies  the  use  of  a  turbulent  slip  velocity  as  a  boundary  condition  for 
the  mean  flow,  the  no-slip  boundary  condition  is  generally  used  except 
for  extremely  thin  boundary  layers. 

At  a  symmetry  boundary  the  normal  gradients  of  k  and  e  are  set  to 
zero.  For  periodic  boundaries,  the  periodicity  condition  is  enforced 
implicitly  using  a  periodic  matrix  solver. 

4.3.2  Initial  Conditions 

Two  methods  can  be  employed  to  generate  initial  conditions  for  the 
k-e  equations.  First,  when  experimental  data  is  available  for  the 
turbulence  intensities,  the  turbulence  kinetic  energy  is  determined 


from: 
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k  =  1/2(  u2  +  v2  +  w2  ) 


(4.15) 


With  the  Reynolds  shear  stress  measured,  the  definition  of  the  eddy 
viscosity  can  be  used  along  with  the  previously  computed  turbulence 
kinetic  energy  to  give  the  turbulence  dissipation  in  the  following  way 
(two  dimensions  are  used  for  example): 

u  v  =  -  Cy  k2/e  Re  3^  Q 

then 

e  =  -Cy  k2  Re  3nQ  /  u  v 

(4.16) 


Thus,  k  and  e  are  known  at  the  inlet  plane  and  the  equations  can  be 
forward-marched  behind  the  mean  flow  equations. 

If  the  experimental  data  is  not  available,  the  algebraic  eddy 
viscosity  model  is  used  to  give  the  eddy  viscosity.  Then  equilibrium  is 
assumed  so  that  production  equals  dissipation  or: 


e  =  Pfc/Re 


(4.17) 


and 

k  =  /  Vpe/(CpRe) 


(4.18) 


A  major  problem  associated  with  the  k-e  model  is  that  the  constants 
in  the  equations  are  based  on  simple  two-dimensional  shear  layers.  In 
addition,  the  turbulence  length  scale  is  generally  over  predicted  which 
artificailly  suppresses  boundary  layer  detachment  in  the  computation. 
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Another  problem  with  the  model  is  that  it  assumes  the  turbulence  is 
isotropic.  True  anisotropy  can  only  be  introduced  by  solving  the  full 
Reynolds  stress  equations  and  not  use  the  eddy  viscosity  hypothesis. 
This  is  computationally  inefficient  especially  for  complex  geometries. 
It  is  unfortunate  that  in  the  field  today,  the  most  simple  models  are 
used  on  the  most  complicated  flows  and  the  most  complex  models  are  used 
on  the  simplest  of  flows.  An  efficient  way  of  introducing  anisotropy  in 
an  approximate  fashion  is  to  assume  that  the  transport  of  the  Reynolds 
stresses  is  very  similar  to  the  transport  of  the  turbulence  kinetic 
energy.  In  this  way,  algebraic  relations  can  be  written  for  the 
Reynolds  stresses. 

4_jl  Algebraic  Reynolds  Stress  Model  ( ARSM) 

For  the  turbomachinery  type  flows  under  consideration  here,  the 
major  generator  of  anisotropy  is  rotation  of  the  flow.  The  rotation 
effects  the  various  Reynolds  stresses  in  an  unequal  way.  The  Coriolis 
force  arising  from  using  a  rotating  coordinate  system  stabilizes  or 
destabilizes  the  turbulence  in  shear  layers.  This  cannot  be  captured  by 
the  standard  k-e  model.  By  relating  the  transport  of  Reynolds  stresses 
to  the  transport  of  k,  Rodi  (1976)  developed  the  following  algebraic 
relation  for  non-rotating  systems. 

u'iu'j=k[2/35ij+(Pij-2P5iJ/3)  ( 1-C2)/(P+e(Cn-1 ) ) ) 

(4.19) 

where  is  the  Kronecker  delta,  P^j  is  the  production  tensor  and 
P=1/3  Pii-  Also,  Ci  =  1.8  and  C2  =  0.6.  Along  similar  lines,  Warfield 
and  Lakshminarayana  (1986)  developed  an  algebraic  relation  for  rotating 
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flows . 


u,iu'j=k[2/3Sij+{(Rij(2-C2)/2^  Pi j-2P«i j/3) ( 1-C2) } /(P+e(CH ) )  ] 

(4.20) 


where 


RiJ=-2flp[e Ju  k  +  +£jpku  iu  k] 


(4.21) 


and  flp  is  the  rotation  vector  and  ejjk  is  the  permutation  tensor. 

Warfield  and  Lakshminarayana  manipulated  equation  4.21  to  give  the 
variable  Cp  Prandtl-Kolmogorov  form  found  in  equation  4.22. 

CM  =  -(2/3)(C2-1)(C2Pk/e+C1-1)/(D1  +D2) 

(4.22) 


where 

Di  =  {(Pk/e)  -  (Ct-1)}2 
D2  =  [4R1(2-C2)/2]2  +  4(C2-1)(2-C2)R2R12 
The  following  natural  groups  appeared  in  their  analysis: 

(Pk/e)  R1=(kfl/e)  R2=3nQ  ^ 

Here,  R2  looks  like  a  gradient  Richardson  number.  Equation  4.22  for 
equilibrium  conditions  with  zero  rotation  gives  a  value  of  Cp  other 
than  the  standard  0.09  thus  it  is  multiplied  by  a  scale  factor  to 
recover  the  standard  value  in  such  conditions.  Also,  when  the  ARSM  is 
applied,  upper  and  lower  limits  are  placed  on  Cp  of  0.2  and  0.025, 
respectively. 
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The  variable  Cy  method  for  applying  the  ARSM  is  stable  and  very 
easy  to  apply  to  existing  k-e  models  and  is  capable  of  modelling  the 
physical  nature  of  rotation  effects  on  turbulence  according  to  Warfield 
and  Lakshminarayana  (1986). 

4.5  Methods  of  Solution 

Integration  of  the  above  turbulence  models  into  the  mean  flow 
equations  is  straightforward.  With  the  solution  to  the  mean  flow 
equations  known  at  station  i,  the  eddy  viscosity  is  determined  from  one 
of  the  various  models.  This  new  eddy  viscosity  is  then  incorporated 
implicitly  in  the  mean  flow  equations  at  station  i+1.  Therefore,  the 
eddy  viscosity  is  lagged  one  streamwise  station  behind  the  mean  flow 
equations.  The  k-e  equations  are  solved  in  an  uncoupled  form  using  ADI 
and  are  marched  in  a  fashion  similar  to  the  mean  flow  equations.  If  the 
ARSM  is  used,  the  k-e  equations  are  computed  at  station  i  using  the 
value  of  Cy  determined  at  i-1.  With  k  and  e  now  known,  a  new  value  of 
Cy  is  determined  from  the  ARSM  and  used  to  find  the  new  value  of  the 
eddy  viscosity.  For  applications  where  the  eddy  viscosity  changes 
drastically  from  one  streamwise  station  to  the  next,  an  iterative  scheme 
may  be  adopted,  however,  this  will  prolong  the  computation  time. 
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CHAPTER  5 

COMPUTATION  OF  COMPLEX  TURBOMACHINERY  FLOWS 

5.1  Introduction 

In  order  to  test  the  prediction  accuracy  of  the  new  parabolic¬ 
marching  method  (NPM),  the  flow  in  several  different  geometries  are 
computed  to  test  the  NPM  method's  ability  to  compute  various  flow 
phenomena.  Initially,  the  turbulent  flow  in  an  S-shaped  duct  is 
computed  since  it  is  characterized  by  a  strong  variation  in  the 
transverse  pressure  gradient.  Next,  the  turbulent  flow  in  an  end  wall 
cascade  is  computed  to  determine  how  well  the  pressure  distribution  and 
end  wall  boundary  layer  growth  is  captured  as  well  as  the  passage 
vortex.  Since  the  wakes  of  turbomachinery  blades  are  very  important  to 
the  designer,  the  turbulent  wake  of  a  cascade  of  airfoils  is  computed. 
This  is  also  a  good  validation  of  the  various  turbulence  models 
incorporated  into  the  NPM  code  since  the  wake  spreading  is  a  strong 
function  of  the  turbulence.  Finally,  the  flow  in  an  axial  flow 
compressor  rotor  is  computed  to  study  the  combined  effects  of  the  blade 
and  end  wall  boundary  layers,  strong  secondary  flow  and  rotation  on  the 
flow  field. 


5.2  Flow  in  an  S-Shaped  Duct 


The  turbulent  flow  in  an  S-shaped  duct  was  computed  using  the  NPM 
method  and  results  compared  to  the  LDV  data  acquired  by  Taylor  et  al. 
(1982).  The  geometry  is  shown  in  Figure  2.9.  No  assumed  pressure  was 
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used  as  initial  conditions  and  the  inlet  flow  was  determined  in  the 
manner  described  in  section  3.8.2.  The  computaional  domain  and  grid 
used  here  is  the  same  as  for  the  laminar  case  presented  in  section 
3.8.2.  The  flow  Reynolds  number  based  on  the  bulk  inlet  velocity  and 
the  duct  width  was  40,000.  This  Reynolds  number  proved  to  be  too  low 
for  successful  application  of  the  k-e  turbulence  model  described  in 
chapter  IV,  thus  the  algebraic  eddy  viscosity  model  was  used. 

For  this  case,  initial  passes  of  the  domain  required  values  of  a 
larger  than  for  the  laminar  flow.  This  trend  was  observed  for  all 
turbulent  flow  computations.  The  convergence  history  for  a=30.0  is 
shown  in  Figure  5.1.  The  chosen  value  of  a  was  the  minimum  value  with 
which  a  convergent  solution  was  obtained.  Typically,  the  more  complex  a 
flow,  the  larger  the  value  of  a.  Nearly  two  orders  of  magnitude  were 
achieved  after  130  global  iterations.  CPU  time  on  the  IBM  3090-200  was 
0.00201  seconds  per  point  per  iteration.  The  convergence  was  not  as 
strong  as  with  the  laminar  flow,  however,  the  same  characteristics  are 
present,  i.e.,  fastest  convergence  for  initial  sweeps  as  small 
wavelength  errors  were  damped  out. 

At  the  initial  measurement  location,  the  computed  streamwise 
velocity  profiles  compare  very  well  with  the  experimental  data, 
especially  near  the  side  walls  (see  Figure  5.2).  Near  the  end  walls, 
the  velocity  gradients  do  not  match  the  data  which  for  internal  flows 
portends  an  incorrect  computed  bulk  pressure  drop.  In  Figure  5.3.  the 
computed  secondary  velocities  at  station  1  is  compared  to  the 
experimental  data  with  very  good  agreement.  Thus  the  inlet  conditions, 
so  important  for  parabolic-marching  methods,  is  accurately  captured. 

At  station  2,  the  computed  streamwise  flow  does  not  compare  well  to 
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the  experimental  data  near  the  center  of  the  duct  in  Figure  5.4. 
Inspection  of  the  profiles  leads  one  to  believe  that  there  may  be  some 
computed  mass  surplus,  however,  the  profiles  very  near  the  side  walls 
were  not  measured  and  so  a  conclusion  regarding  this  cannot  be  drawn. 
The  computed  mass  flow  did  remain  constant  at  each  station  at 
convergence  with  less  than  a  3J  variation. 

The  secondary  flow  at  station  2  has  developed  properly  as  seen  in 
Figure  5.5.  The  inviscid  secondary  flow  is  slightly  overpredicted  yet 
the  maximum  flow  turning  in  the  end  wall  viscous  region  is  well 
predicted.  At  station  4  in  the  second  bend  of  the  duct,  see  Figure  5.6, 
the  computed  screemwise  velocity  profiles  compare  as  well  to  the  data  as 
at  station  2.  At  Y/D=0.1,  however,  the  computed  profile  does  not  match 
the  data  qualitatively  but  at  Y/D=0.9  there  is  excellent  agreement 
between  the  computation  and  the  experimental  data  of  Taylor  et  al. 
(1982).  The  computed  secondary  flow  at  this  station,  shown  in  Figure 
5.7,  is  in  excellent  agreement  with  the  data  at  all  measurement 
locations.  Finally  in  Figure  5.8,  the  computed  streamwise  velocity 
profiles  at  station  5  continue  to  show  only  adequate  agreement  with  the 
data  near  the  center  of  the  duct.  The  best  predictions  come  near  the 
side  walls.  In  Figure  5.9,  the  computed  secondary  flow  profiles  compare 
well  to  the  experimental  data  except  near  the  side  walls  especially  at 
Y/D=0.9.  Here,  the  measured  transverse  velocity  has  an  s-shaped  profile 
and  the  magnitude  of  the  velocity  near  the  end  wall  is  large.  The 
computed  velocity  profile  does  not  exhibit  this  behavior.  This  strong 
secondary  flow  may  have  been  missed  due  to  the  coarse  nature  of  the 
computational  grid.  Again  it  should  be  noted  that  the  velocity 
gradients  very  near  the  wall  do  not  compare  well  to  the  data  which  may 
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give  rise  to  a  poor  prediction  of  the  bulk  pressure  drop.  Indeed,  the 
computed  transverse  pressure  gradient  is  well  predicted  everywhere  while 
the  bulk  pressure  drop  is  not  well  predict  as  seen  in  Figure  5.10.  With 
the  excellent  agreement  with  the  data  for  the  laminar  case  for  this  same 
geometry,  the  less  than  stellar  prediction  of  the  turbulent  flow 
implicates  the  turbulence  model.  The  value  of  the  near  wall  velocity 
gradient  is  strongly  dependent  upon  the  turbulence  model.  Good 
predictions  of  the  skin  friction  in  two  dimensions  have  been  achieved  by 
many  users  of  the  Baldwin  and  Lomax  (1978)  algebraic  model.  Therefore, 
simple  extension  of  the  two-dimensional  model  to  three  dimensions 
without  regard  for  the  actual  physics  of  a  three-dimensional  boundary 
layer  may  be  the  reason  for  the  reduced  accuracy.  The  next  challenge 
then  is  to  compute  the  turbulent  flow  in  a  turbomachinery  blade  row. 

5.3  Flow  in  an  End  Wall  Cascade 

The  three-dimensional  turbulent  flow  in  a  turbine  end  wall  cascade 

was  computed  using  the  NPM  method.  The  flow  was  measured  using  a  hot 

wire  probe  by  Flot  and  Papailiou  (1975).  The  cascade  was  made  up  of 

o 

NACA  65-12-Aiq-IO  blades  with  a  stagger  angle  of  -15  .  The  span  to 
chord  ratio  was  2.1,  the  pitch  to  chord  ratio  was  0.8  and  the  flow 
Reynolds  number  was  389,000.  The  geometry  and  measurement  locations  are 
given  in  Figure  5.11.  Only  the  blade  region  was  computed  and  the 
computational  grid  consisted  of  30  points  in  the  streamwise  direction 
and  a  23x23  grid  in  the  cross  plane.  Both  the  algebraic  eddy  viscosity 
model  and  the  k-e  turbulence  model  were  used  to  compute  this  flow. 
There  was  very  little  difference  in  computational  results  between  the 
two  models  but  the  use  of  the  k-e  model  required  more  computation  time. 
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The  following  results  were  from  the  two-equation  model.  The  no-slip 
boundary  condition  was  used  on  the  blade  surfaces  as  well  as  on  the  end 
wall  boundary.  A  symmetry  boundary  condition  was  applied  at  the  mid¬ 
span.  The  inviscid  velocity  distribution  was  used  as  inlet  conditions 
at  the  blade  leading  edge  to  start  the  computation. 

The  convergence  history  is  presented  in  Figure  5.12.  For  this  case, 
a  was  set  to  20  and  the  pressure  residuals  dropped  one  and  one  half 
orders  of  magnitude  in  50  global  iterations.  At  this  point,  the 
convergence  flattened  out  and  the  computation  was  stopped  in  order  to 
conserve  CPU  time. 

The  computed  pressure  distribution  at  the  mid-span  is  compared  to 
the  experimental  data  in  Figure  5.13  with  excellent  agreement.  The 
inviscid  pressure  distribution  computed  from  the  panel  code  of  Giesing 
(1964),  used  as  the  initial  pressure  distribution,  is  also  presented  in 
Figure  5.13.  The  pressure  near  the  end  wall  is  very  similar  to  that  at 
the  mid-span  which  gives  creedence  to  the  assumption  that  the  normal 
gradient  of  pressure  is  zero  following  classical  boundary  layer  theory. 
The  good  agreement  with  the  experimental  data  indicates  that  the  NPM 
method  is  capable  of  computing  the  viscous  pressure  field  accurately 
unlike  the  PL  method  and  various  other  space-marching  methods. 

The  computed  streamwise  end  wall  velocity  profiles  at  44?  chord  and 
88?  chord  are  compared  to  the  experimental  data  in  Figures  5.14  and 
5.15,  respectively.  The  agreement  with  the  data  is  very  good  near  the 
suction  side  and  near  mid  pitch  at  both  streamwise  locations.  Near  the 
pressure  side,  the  results  are  not  as  good  and  this  may  be  due  to  a  poor 
prediction  of  the  pressure  side  blade  boundary  layer.  Unfortunately,  no 
blade  boundary  layers  were  measured  to  verify  this.  Figure  5.16  shows 
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the  computed  streamwise  momentum  thickness  compared  to  the  the 
experimental  data.  The  blade  to  blade  variation  is  well  predicted  as  is 
the  level  of  the  momentum  thickness.  One  can  see  that  the  secondary 
flow  is  convecting  the  low  momentum  fluid  to  the  suction  side  from  the 
pressure  side  boundary  layers  thereby  thinning  them.  The  computed 
secondary  end  wall  velocity  profiles  are  compared  to  the  experimental 
data  at  44%  chord  and  88?  chord  in  Figures  5.17  and  5.18,  respectively. 
At  all  locations,  the  overturning  of  the  flow  near  the  end  wall  is 
reasonably  well  predicted  while  the  flow  underturning  in  the  outer 
reaches  of  the  boundary  layer  is  underpredicted.  Since  the  streamwise 
profiles  were  well  predicted,  one  can  assert  that  the  error  is  due  to  a 
poor  prediction  of  the  transverse  pressure  gradient  near  the  end  wall 
which  drives  the  secondary  flow.  Overall,  however,  the  prediction  of 
the  secondary  flow  is  good. 

The  passage  vortex  is  clearly  visable  from  a  plot  of  the  computed 
secondary  velocity  vectors  in  Figures  5.19  through  5.21.  Careful 
inspection  reveals  that  the  center  of  the  passage  vortex  is  convected 
from  the  pressure  side  to  the  suction  side  boundary.  Finally,  the 
computed  mass  averaged  secondary  flow  kinetic  energy  ksec  is  presented 
in  Figure  5.22.  The  secondary  flow  kinetic  energy  may  be  thought  of  as 
unrecoverable  energy  thus  it  is  indicative  of  the  viscous  flow  losses 
within  the  passage.  From  the  results  presented  in  this  section,  it  is 
apparent  that  the  NPM  method  is  capable  of  computing  the  viscous  flow  in 
the  end  wall  region  of  a  turbomachinery  blade  row.  The  major  feature  of 
such  a  flow,  namely  the  passage  vortex  system,  is  well  predicted. 
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Figure  5.20  Secondary  Flow  Velocities  at  66%  Chord  for  the 
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Figure  5.21  Secondary  Flow  Velocities  at  88*  Chord  for  the 
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Figure  5.22  Mass  Averaged  Secondary  Flow  Kinetic  Energj 


for  the  End  Wall  Cascade 


152 


5.4  Compressor  Cascade  Wake  Flow 

The  new  parabolic-marching  method  was  used  to  compute  the  turbulent 

wake  of  a  compressor  cascade.  The  flow  field  was  measured  by  Hobbs  et 

al.  (1980)  using  a  hot  film  probe  for  the  near  wake  and  a  five  hole 

probe  for  the  far  wake.  The  blades  are  double  circular  arc  airfoils 

0 

with  a  stagger  angle  of  20.77  .  The  space  to  chord  ratio  is  0.6  and 
the  flow  Reynolds  number  is  588,000.  Although  the  experimental  facility 
used  end  wall  suction  to  ensure  two-dimensionality,  the  computation  was 
carried  out  assuming  an  end  wall  was  present  and  only  results  at  the 
mid-span  were  considered.  The  geometry  and  computational  grid  at  the 
mid-span  is  shown  in  Figure  5.23.  The  inviscid  pressure  computed  from 
the  panel  method  code  of  Giesing  (1964)  was  used  as  the  initial  pressure 
field.  The  corresponding  inviscid  velocity  at  the  inlet  was  used  as 
inlet  conditions  as  was  done  for  the  Plot  and  Papailiou  (1975)  end  wall 
cascade  in  section  5.3.  In  order  to  save  CPU  time  and  since  the  blade 
boundary  layers  are  very  thin,  the  inviscid  inlet  condition  was  applied 
near  the  mid-chord  region  where  the  boundary  layer  thickness  of  a 
typical  turbulent  boundary  layer  reached  a  height  equal  to  the  distance 
of  the  first  grid  point  from  the  surface.  The  no-slip  condition  was 
used  on  the  blade  surfaces  and  the  slip  condition  was  applied  on  the  end 
wall  in  an  effort  to  achieve  two-dimensionality. 

The  flow  was  computed  with  both  the  algebraic  eddy  viscosity  model 
and  the  two-equation  turbulence  model.  When  the  algebraic  model  was 
used,  the  best  convergence  was  achieved  with  a  equal  to  10  as  seen  in 
Figure  5.24.  The  wake  centerline  velocity  computed  with  the  algebraic 
eddy  viscosity  model  is  compared  to  the  experimental  data  in  Figure 
5.25.  x  is  the  distance  from  the  trailing  edge  and  Bx  is  the  axial 
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chord  length.  Also  displayed  is  the  wake  correlation  of  Raj  and 
Lakshrainarayana  (1973).  The  agreement  between  the  computation, 
experimental  data,  and  correlation  is  good,  however,  the  computed  wake 
defect  is  underpredicted  in  the  near  wake  region  while  it  is 
overpredicted  in  the  far  wake  region.  The  comparison  of  the  computed 
wake  profiles  with  the  experimental  data  is  given  in  Figure  5.26.  In 
the  very  near  wake,  the  computed  velocity  profiles  do  not  compare  well 
with  the  experimental  data  near  the  wake  centerline.  In  this  region, 
the  velocity  gradients  are  very  large  and  the  NPM  method  had  difficulty 
computing  them.  In  the  outer  flow  region,  the  velocity  profile  is  very 
well  predicted  owing  to  the  strong  influence  of  the  boundary  layer  which 
was  reasonably  well  captured.  The  nature  of  the  turbulence  is  very 
complicated  near  the  centerline  of  the  wake  and  the  algebraic  turbulence 
model  does  not  seem  to  have  the  power  to  resolve  the  strong  interaction. 

The  far  wake  is  only  adequately  predicted.  The  wake  centerline 
position  is  below  the  measured  position  and  the  wake  spreading  is  less 
than  measured.  The  wake  spreading  is  a  strong  function  of  the 
turbulence  and  the  ad  hoc  approximations  made  in  the  algebraic 
turbulence  model  may  contribute  to  the  poor  prediction. 

The  two-equation  k-e  model,  which  includes  more  physical  properties 
in  its  formulation,  was  used  in  an  effort  to  improve  the  prediction. 
The  convergence  history  shown  in  Figure  5.27  is  almost  identical  to  that 
using  the  algebraic  eddy  viscosity  model.  The  extra  physics  included  in 
the  k-e  model  have  indeed  improved  the  prediction  of  the  wake  flow. 
The  computed  centerline  velocity,  shown  in  Figure  5.28,  indicates  the 
improved  prediction  accuracy  especially  in  the  far  wake  region.  The 
wake  defect  is  again  underpredicted  in  the  near  wake  region  by  a 
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substantial  amount,  however,  in  the  far  wake,  the  defect  matches  the 
experimental  data  almost  exactly.  The  spreading  rate  is  also  accurately 
predicted  in  the  far  wake. 

The  wake  velocity  profiles  computed  with  the  k-e  model  are  shown  in 
Figure  5.29.  Again  the  boundary  layer  profiles  are  reasonably  well 
predicted,  however,  the  near  wake  profiles  show  the  same  lack  of 
agreement  with  the  centerline  data  as  with  the  computation  with  the 
algebraic  model.  In  the  far  wake,  the  advantages  of  the  k-e  model  are 
clear  with  excellent  agreement  between  the  computation  and  the 
measurements.  The  wake  centerline  position  as  well  as  the  wake 
spreading  is  accurately  predicted  even  though  the  near  wake  predictions 
were  not  as  accurate.  This  leads  one  to  believe  that  the  far  wake  is 
not  strongly  influenced  by  the  near  wake.  Overall,  the  prediction  using 
the  k-e  model  is  superior  to  the  prediction  using  the  algebraic  model 
for  all  regions  of  the  flow.  The  computed  turbulence  intensities  are 
compared  to  the  experimental  data  for  the  near  wake  region  in  Figure 
5.30.  On  the  the  blade  surface,  the  turbulence  intensities  compare  very 
well  with  the  data.  This  is  not  surprising  since  the  established  k-e 
model  has  been  'fine  tuned'  for  wall  bounded  shear  layers.  In  the  near 
wake  region,  the  peak  intensities  are  not  well  captured.  This  can  be 
attributed  to  the  overprediction  of  the  turbulence  dissipation.  For 
some  coarse  grid  computations,  this  turbulence  dissipation  was  so  large 
that  no  turbulence  could  be  sustained  past  the  near  wake  region.  Thus, 
there  appears  to  be  some  grid  dependence  for  this  wake  flow.  It  should 
be  noted,  though,  that  for  a  successful  wake  computation,  special  care 
had  to  be  taken  to  ensure  that  the  flow  did  not  separate  at  the  trailing 
edge.  Such  a  situation  generated  incorrect  turbulence  kinetic  energy 


and  turbulence  dissipation  distributions  in  the  strong  interaction 
region  which  generally  leads  to  divergent  behavior. 

The  computed  pressure  distribution  in  the  wake  region  is  shown  in 
Figure  5.31.  The  most  important  observation  is  that,  while  the 
streamwise  pressure  gradient  in  the  outer  flow  is  nearly  constant,  the 
pressure  drops  dramatically  near  the  wake  centerline  close  to  the 
trailing  edge.  This  localized  effect  is  due  to  the  acceleration  of  the 
wake  centerline  velocities.  Thus,  the  transverse  pressure  gradient  is 
not  constant  across  the  shear  layer  in  the  near  wake  region  even  for  a 
non  curving  wake. 

5.5  Flow  in  an  Axial  Flow  Compressor  Rotor 

The  final  test  of  the  NPM  method  was  on  the  PSU  axial  flow 
compressor  rotor.  The  facility  is  described  in  depth  by  Lakshminarayana 
(1980).  The  hub  to  tip  ratio  is  0.5  with  the  radius  of  the  annulus  wall 

being  0.466  meters.  The  rotor  is  made  up  of  21  NACA  65-010  blades  with 

0  0 

the  stagger  angle  varying  from  22.5  at  the  hub  to  45  at  the  tip.  The 
overall  performace  of  this  compressor  is  given  in  Figure  5.32.  The 
computations  performed  here  were  for  the  design  flow  coefficient  of 
4>=0.56.  Although  the  tip  clearance  region  is  roughly  0.15  cm,  the  tip 
clearrnce  effects  were  not  included  in  the  computation.  Also,  the 
measured  relative  inlet  Mach  number  is  0.085  and  the  Mach  number  based 
on  tip  speed  is  0.153  thus  the  assumption  of  incompressibility  is  valid. 

The  inviscid  inlet  conditions  at  the  leading  edge  and  the  initial 


inviscid  pressure  field  were  generated  by  stacking  several  two- 
dimensional  panel  solutions  using  Oiesing's  (1964)  program.  Turbulent 
slip  boundary  conditions  were  used  on  all  blade  surfaces.  Major 
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Figure  5.31  Pressure  Distribution  for  Hobbs  Cascade,  Alg.  Model 
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difficulties  were  encountered  when  the  k-e  turbulence  model  was  applied 
to  this  case,  thus  the  computation  was  performed  using  the  algebraic 
eddy  viscosity  model. 

The  convergence  history  for  the  rotor  computation  using  a  equal  to 
25  can  be  seen  in  Figure  5.33.  Convergence  of  one  and  one  half  orders 
of  magnitude  was  achieved  in  100  global  iterations.  The  computation  was 
stopped  to  conserve  CPU  time.  The  computed  pressure  distribution  at 
various  radial  locations  is  compared  to  the  experimental  data  of  Sitaram 
and  Lakshminarayana  (1983)  in  Figures  5.34  through  5.37.  The  suction 
peak  is  well  captured  near  the  mid-span,  however,  it  is  overpredicted  in 
the  hub  wall  region.  The  suction  side  pressure  distribution  agrees 
favorably  with  the  experimental  data.  The  pressure  side  Cp  agrees  well 
with  the  data  at  a  hub  to  tip  ratio  R=0.832,  however,  at  lower  values  of 
R,  a  hump  in  the  computed  distribution  is  present  near  the  trailing 
edge.  This  may  have  arisen  from  the  fact  that  the  wake  region  has  not 
been  computed  and  without  the  interaction  of  the  wake,  the  flow  turns 
very  sharply  (locally)  to  follow  the  blade  trailing  edge  rather  than 
following  the  outer  inviscid  flow.  Such  turning  tends  to  increase  lift 
which  is  evident  in  the  computed  pressure  distribution. 

The  computed  suction  side  boundary  layer  profiles  are  compared  to 
the  hot  wire  data  of  Pouagare  et  al.  (1985)  in  Figure  5.38.  The 
agreement  with  the  data  is  good  near  the  hub  and  mid-span  regions.  At 
R=0.918,  the  prediction  accuracy  is  less  than  adequate.  It  is  important 
to  remember  that  the  tip  clearance  effect  has  not  been  included  in  the 
computation  and  the  measured  data  at  this  location  may  indeed  include 
the  interaction  of  the  leakage  Jet  with  the  blade  boundary  layer.  Near 
the  trailing  edg*»,  the  boundary  layer  profiles  are  not  well  predicted 
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Figure  5.35  Pressure  Distribution  at  R=.670  for  the  Rotor 


Figure  5.36  Pressure  Distribution  at  R=.750  for  the  Rotor 
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which  may  be  due  to  the  lack  of  wake  interaction  in  the  computation. 

The  computed  pressure  side  boundary  layer  profiles  are  compared  to 
the  experimental  data  in  Figure  5.39.  Here,  the  agreement  is  not  good. 
These  boundary  layers  are  so  thin  that  the  measurements  must  be  suspect 
in  the  error  analysis.  Even  so,  the  computed  profiles  do  not  show  the 
proper  character  at  R=0.750.  This  may  be  due  to  the  poor  pressure 
prediction.  The  computed  radial  velocity  profiles  at  S=.99  for  the 
pressure  and  suction  side  are  compared  to  the  data  in  Figure  5.40.  The 
results  are  qualitatively  correct,  however,  the  magnitudes  are  somewhat 
underpredicted.  This  may  be  due  to  the  the  use  of  the  turbulent  slip 
boundary  condition.  The  radial  flow  near  the  blade  surface  is  driven  by 
the  normal  gradient  of  the  streamwise  velocity.  Thus,  since  the  slip 
condition  effectively  reduces  this  gradient,  the  radial  flows  will  not 
develop  as  is  the  case.  The  computed  hub  wall  boundary  layer  is 
compared  to  the  experimental  data  in  Figure  5.41.  The  boundary  layer 
growth  on  the  hub  wall  is  well  captured  as  is  the  actual  profile.  The 
computed  secondary  flow  vectors  are  compared  to  the  experimental  data  of 
Mur thy  and  Lakshminarayana  (1987)  in  the  hub  wall  region  in  Figures  5.42 
through  5.44.  Overall,  the  secondary  flow  has  the  correct  trend  near 
the  hub  wall,  however,  the  magnitude  is  slightly  overpredicted. 
Although  artificial,  due  to  the  lack  of  tip  clearance  modelling,  a 
passage  vortex  can  be  seen  developing  near  the  annulus  wall  in  Figure 
5.42a.  Radial  flows  are  small  at  this  streamwise  station.  At  S=.58,  in 
Figure  5.43a,  the  radial  flows  are  stronger  as  are  the  secondary  flows. 
The  passage  vortex  near  the  annulus  wall  is  more  strongly  developed. 
Finally,  the  radial  flows  are  so  strong  at  S=1.0,  see  Figure  5.44a,  that 
much  of  the  mid-span  momentum  is  convected  to  the  blade  tip.  This  may 
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)  Figure  5.42  Secondary  Flow  Vectors  at  X=.35for  the  Rotor 
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b)  Measured:  Mur  thy  and  Lakshminarayana  (1987) 
i  Figure  5.42  Secondary  Flow  Vectors  at  X=.35  for  the  Rotor 
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Figure  5.43  Secondary  Flow  Vectors  at  X=.58  for  the  Rotor 
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b)  Measured:  Murthy  and  Lakshminarayana  (1987) 

Figure  5.43  Secondary  Flow  Vectors  at  X=.58  for  the  Rotor 
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b)  Measured:  Murthy  and  Lakshminarayana  (1987) 

Figure  5. 44  Secondary  Flow  Vectors  at  X=1.0  for  the  Rotor 
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be  another  reason  for  the  poor  boundary  layer  predictions  at  S=.99.  The 
interaction  of  the  tip  leakage  Jet  tends  to  slightly  reduce  these  radial 
flows  so  perhaps  a  more  accurate  computation  can  be  achieved  with  some 
tip  clearance  modelling.  In  the  hub  wall  region,  however,  the 
comparison  with  the  data  is  very  good  at  this  streamwise  station. 
Overall,  the  major  features  of  the  rotor  flow  have  been  captured 
accurately. 

Finally,  the  computed  stagnation  pressure  loss  coefficient,  4>ioss» 
based  on  the  blade  tip  speed  is  compared  to  the  experimental  data  at  the 
exit  in  Figure  5.45.  The  agreement  with  the  data  is  surprisingly  good. 
The  losses  on  the  pressure  side  are  somewhat  underpredicted  and  the 
computed  loss  core,  near  n/c=0.3,  is  closer  to  the  suction  side  than  the 
measured  location  of  n/c=0.4.  Still,  the  stagnation  pressure  loss  is 
one  of  the  most  difficult  flow  quantities  to  capture  and  the  NPM  method 
has  done  a  reasonably  good  job  of  computing  it. 
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CHAPTER  6 


CONCLUDING  REMARKS 


6.1  Summary  of  Conclusions 


Existing  single  pass  space-marching  methods  for  solving  the 
partially  parabolized  Navier-Stokes  equation  have  various  drawbacks 
associated  with  their  use  on  complex  flows.  Several  require  Poisson 
correction  equations  which  are  cumbersome  to  use  and/or  do  not  yield  the 
viscous  pressure  field.  One  method  in  particular  due  to  Pouagare  and 
Lakshminarayana  (1986),  refered  to  as  the  PL  method,  overcomes  many  of 
these  drawbacks.  The  PL  method,  however,  was  found  through  analysis  to 
be  unconditionally  unstable  when  used  as  an  iterative,  multiple  pass, 
parabolic-marching  method.  This  restricts  its  usefullness  in  predicting 
flows  with  strong  pressure  interaction.  The  PL  method  was  capable  of 
computing  many  mildly  elliptic  internal  flows  because  the  continuity 
equation  was  directly  coupled  to  the  momentum  equations  in  the  solution 
process  and  the  streamwise  pressure  gradient  condition  was  relaxed. 
This  relaxation  of  the  streamwise  pressure  gradient  was  the  reason  for 
the  global  instability,  where  global  refers  to  multiple  passes  of  the 
domain. 

A  modification  was  made  to  the  PL  method  in  an  effort  to  improve  its 
global  stability  characteristics.  The  modified  PL  method,  refered  to  as 
the  MPL  method,  relaxed  the  transverse  pressure  gradients  still  keeping 
the  continuity  equation  coupled  to  the  momentum  equations.  Analysis 
showed  that  this  system  was  stable  for  a  forward-marching  integration  of 


the  equation  yet  globally  unstable  for  multiple  passes  of  the  domain. 
Clearly,  for  global  stability,  the  pressure  gradients  must  not  be 
relaxed.  Since  some  condition  must  be  relaxed  for  the  solution  of  the 
coupled  form  of  the  Navier-Stokes  equation  on  complex  geometries,  the 
continuity  equation  mas  determined  to  Qe  the  appropriate  condition  to 
relax. 

Following  pseudocompressibility  theory,  a  pressure  term  was 
introduced  into  the  continuity  equation  giving  a  system  which  appears 
similar  to  Chorin's  (''967)  method  but  with  no  timewise  gradients.  The 
new  parabolic-marchii  g  method  or  NPM  method  was  found  to  be  stable  when 
a  forward-marching  integration  was  used  used  as  a  solution  procedure. 
Most  importantly,  the  NPM  method  was  found  to  be  unconditionally  stable 
for  a  multiple  pass,  global  iteration  procedure.  The  value  of  the 
relaxation  parameter  in  the  continuity  equation  was  found  to  vary  from 
zero  to  thirty  with  a  general  trend  of  increasing  value  for  increasing 
flow  complexity.  In  addition,  the  mass  flow  was  found  to  be  conserved 
at  convergence,  nowever,  convergence  was  generally  slow  as  is  typical  of 
parab~lic-marching  methods. 

The  LBI  scheme  used  to  solve  the  equation  in  the  cross  plane  was 
found  to  introdu'  an  inconsistent  splitting  error  when  used  in 
parabolic-marching  methods.  This  splitting  error  is  strongly  dependent 
upon  the  way  in  which  the  main  vector  operator  in  the  LBI  scheme  is 
split  and  reduces  accuracy  significantly.  Depending  upon  the  flow  under 
consideration,  use  of  one  splitting  order  over  the  other  can  lead  to 
divergent  behavior.  This  was  demonstrated  for  the  laminar  flow  in  a 
strongly  curved  duct.  A  new  iterative  solution  procedure  was  .developed 
which  reduces  the  effect  of  the  splitting  error  and  gives  a  result  which 
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maintains  the  formal  second  order  accuracy  of  the  discretization. 

The  NPM  method  using  the  new  solution  procedure  was  calibrated  by 
computing  the  developing  flow  in  a  straight  square  duct  as  well  as  the 
laminar  flow  in  an  S-shaped  square  duet.  The  agreement  with  the 
analysis  and  the  experimental  data  was  excellent  in  all  regions  of  the 
flows.  Since  the  streamwise  diffusion  of  momentum  was  neglected  to  give 
the  partially  parabolized  Navier-Stokes  equation,  the  flow  'sees'  a 
change  in  curvature  through  the  pressure  field  only.  Without  proper 
transmittion  of  pressure  ellipticity,  accurate  comparisons  with  the  data 
are  impossible.  The  results  presented  in  this  work  are  very  accurate 
for  the  laminar  flow;  thus,  the  S-shaped  duct  computation  showed  that 
pressure  ellipticity  could  be  transmitted  effectively  using  the  NPM 
method.  Further  validation  of  the  NPM  method  was  carried  out  on  several 
complex,  turbulent  flows.  The  turbulent  flow  in  an  S-shaped  duct  was 
computed  with  good  agreement  with  the  available  experimental  data.  The 
bulk  pressure  drop  was  overpredicted  but  the  transverse  gradients  due  to 
the  duct  curvature  were  well  captured.  The  error  in  the  bulk  pressure 
drop  may  be  due  to  a  poor  prediction  near  the  walls.  When  the  accuracy 
exhibited  for  the  laminar  flow  is  considered,  one  must  suspect  the 
extension  to  three  dimensions  of  the  two-dimensional  Baldwin  and  Lomax 
(1978)  algebraic  eddy  viscosity  model  as  one  reason  for  the  error. 

The  two  equation  k-e  turbulence  model  was  used  to  compute  the  flow 
in  a  rectilinear  end  wall  turbine  cascade.  The  end  wall  flow 
computation  compared  very  well  to  the  experimental  data.  The 
development  of  the  passage  vortex  was  accurately  captured  as  was  the 
pressure  distribution.  The  wake  of  a  similiar  cascade  was  computed 
using  both  the  algebraic  eddy  viscosity  model  and  the  k-e  model.  The 
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computed  near  wake  flow  compared  very  well  to  the  experimental  data  in 
the  outer  region  and  only  adequately  near  the  wake  centerline.  The  far 
wake,  however,  was  not  well  predicted  with  the  algebraic  eddy  viscosity 
model.  Indeed,  both  the  location  of  the  wake  centerline  and  the  wake 
decay  did  not  match  the  data.  Recall  that  the  wake  spreading  is  a 
strong  function  of  the  turbulence  and  the  degree  of  empericism  and  lack 
of  physical  considerations  found  in  the  algebraic  turbulence  model 
yielded  a  spreading  rate  that  was  too  slow.  On  the  other  hand,  when  the 
k-e  model  was  used,  the  far  wake  computation  was  extremely  accurate. 

Finally,  the  flow  in  an  axial  flow  compressor  rotor  was  computed 
using  the  NPM  method.  This  particular  flow  proved  to  be  very  difficult 
to  compute.  The  algebraic  eddy  viscosiy  model  was  used  for  the  simple 
reason  that  both  the  k-e  model  and  the  ARSM  showed  divergent  behavior 
for  reasons  unknown  to  the  author.  Perhaps  the  forward-marching 
solution  procedure  used  for  the  k  and  e  equations  is  not  appropriate 
for  this  type  of  flow.  The  suction  side  boundary  layers  as  well  as  the 
pressure  distribution  were  well  predicted.  The  computed  pressure  side 
boundary  layers  and  the  exit  radial  flows  did  not  compare  well  with  the 
experimental  data.  The  lack  of  accuracy  for  the  radial  flows  may  be 
attributed  to  the  use  of  the  turbulent  slip  condition  which  tends  to 
reduce  che  near  wall  normal  gradient  of  the  streamwise  velocity.  This 
gradient  is  the  main  mechanism  driving  the  radial  flows  near  the  blade 
surfaces.  The  secondary  flows  were  well  predicted  near  the  hub  wall  as 
were  the  stagnation  pressure  losses  at  the  rotor  exit. 

Overall,  the  NPM  method  showed  its  usefullness  in  computing 
turbulent  turboraachinery  flows.  The  rotor  computation  will  require  more 
modelling,  a  finer  grid,  and  integration  of  the  wake  flow  if  a  more 
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accurate  prediction  is  to  be  achieved. 

6.2  Suggestions  for  Future  Research 

Future  research  into  the  application  of  the  NPM  method  should 
include  the  development  of  a  method  for  accelerating  the  convergence  of 
the  scheme.  Many  iterative  methods  have  incorporated  multigrid  type 
algorithms  to  enhance  convergence.  Such  a  method  may  be  successfully 
incorporated  into  the  NPM  method.  Since  only  pressure  information  needs 
to  be  quckly  propagated  from  downstream  to  upstream,  an  approximate  form 
of  the  Poisson  equation  for  the  pressure,  similar  to  the  one  used  by 
Tenpas  and  Pletcher  (1987),  may  be  sufficient  to  significantly  improve 
the  convergence. 

Another  area  where  the  NPM  method  can  be  improved  is  the  iterative 
LBI  method  employed  at  each  streamwise  station.  There  may  be  faster 
iterative  schemes  available  which  will  yield  the  same  result. 

Finally,  the  area  which  needs  the  most  attention  and  which  is 
vitally  important  to  the  successful  computation  of  complex  flows  is 
turbulence  modelling,  or  rather  the  application  of  existing  turbulence 
models  to  the  unique  problems  in  turbomachinery.  Often,  ad  hoc 
approximations  are  required  in  order  to  achieve  a  converged  solution. 
These  approximations,  e.g.,  location  of  transition,  have  a  substantial 
effect  on  the  outcome  of  a  computation.  What  is  required  is  a 
systematic  study  of  the  established  turbulence  models  and  their 
extension  to  three-dimensional  turboraachinery  flows. 
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APPENDIX 


COMPUTATIONAL  GRID  GENERATION 


The  governing  equation  must  be  solved  on  discrete  points  in  the 
physical  domain.  For  logistical  reasons,  the  distribution  of  points  is 
not  random.  Applying  the  boundary  conditions  is  the  most  complicated 
aspect  of  the  numerical  method  and  is  simplified  by  writing  the 
equations  in  generalized  coordinates.  What  is  now  required  is  to 
generate  the  distribution  of  points  in  the  physical  domain  which  is 
aligned  with  the  generalized  coordinates. 

There  are  three  types  of  point  distributions  or  grids.  They  are  0- 
grids,  C-grids,  and  H-grids.  For  a  parabolic  solution  procedure  for  the 
mean  flow  equations,  the  0-grids  and  C-grids  cannot  be  used.  This  is 
because  the  grid  line  along  which  the  equations  are  marched  is  not 
necessarily  aligned  with  the  streamwise  direction.  For  both  grids,  in 
the  leading  edge  region,  the  marching  direction  is  perpendicular  to  the 
streamwise  direction.  This  is  also  the  case  in  the  trailing  edge  region 
for  0-grids.  Thus,  only  H-grids  can  be  used  for  the  computation  of 
turbomachinery  flows  with  a  parabolic  marching  scheme. 

There  are  several  methods  for  generating  these  H-grids  which  can  be 
grouped  into  partial  differential  equation  methods  and  algebraic 
methods.  The  PDE  method,  of  course,  solves  partial  differential 
equations  to  generate  the  distribution  of  computational  points.  With 
such  methods,  the  user  has  great  control  of  orthogonality  of  the  grid 
lines  in  the  physical  domain.  With  algebraic  methods,  orthogonality  of 
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the  grid  lines  is  not  always  assured.  In  fact,  the  grid  lines  can  be 
highly  sheared  for  periodic  grids  for  cascade  geometries  with  large 
stagger.  This  shearing  of  the  grid  is  detrimental  to  the  accuracy  of 
partially  parabolized  methods  since  the  second  order  cross  derivatives 
with  respect  to  the  streamwise  direction  are  neglected.  These  cross 
derivatives  are  very  important  in  viscous  flows  with  a  non-orthogonal 
grid.  However,  if  a  PDE  method  is  used  to  generate  the  grid,  large 
bulges  appear  in  the  leading  and  trailing  edge  regions  due  to  the 
enforcement  of  orthogonality.  Experience  has  shown  that  this  leads  to 
numerical  problems  which  inhibit  the  computation  of  the  flow. 
Therefore,  an  algedraic  grid  is  used  here  for  the  computation  of 
turbomachinery  flows. 

fhe  following  describes  the  method  used  here  to  generate  three- 
dimensional  grids  for  a  turbomachinery  rotor  in  the  R-0-Z  coordinate 
system.  First,  pseudo-streamsheets  must  be  generated,  i.e.,  two- 
dimensional  grids  in  the  R-0  plane.  These  streamsheets  are  refered  to 
as  'pseudo'  since  they  are  determined  by  the  geometry  of  the  hub  and 
annulus  walls  and  not  by  the  flow.  With  the  hub  boundary  and  annulus 
boundary  defined  from  the  input  data,  an  arbitrary  number  of  points  is 
chosen  along  the  machine  axis  to  define  the  pseudo-streamsheets.  A 
uniform  distribution  is  sufficient.  Lines  are  then  extended  from  these 
points  through  the  hub  terminating  at  the  annulus.  The  corresponding 
values  of  R  and  Z  for  each  line  are  determined  through  bilinear 
interpolation  using  the  input  data. 

With  these  lines  now  generated,  a  distribution  of  points  in  the  Z- 
direction  is  chosen.  Typically,  exponential  stretching  is  used  to  pack 
grid  points  near  the  walls  where  large  flow  gradients  are  anticipated. 
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When  completed  for  each  line,  the  result  is  a  two-dimensional  grid  in 
the  R-Z  plane  where  the  streamwise  grid  lines  are  the  cylindrical 
pseudo-streamsheets.  Two-dimensional  H-grids  are  now  generated  on  the 
pseudo-streamsheets  using  the  following  procedure. 

First,  a  new  blade  profile  is  interpolated  from  the  input  data  for  a 
particular  streamsheet.  This  new  profile  is  then  split  into  a  lower  and 
upper  half  and  the  lower  half  is  shifted  in  the  negative  8  direction  by 
an  amount  equal  to  2n/number  of  blades.  This  provides  an  upper  and 
lower  bound  for  the  grid.  Next,  the  upper  and  lower  boundaries  are 
extended  upstream  and  downstream  of  the  blades  in  the  Z-direction  to  the 
predetermined  extents  of  the  domain.  Now  the  axial  distribution  of 
points  must  be  generated.  Once  again,  exponential  stretching  is  used  to 
pack  grid  points  near  the  leading  and  trailing  edges  of  the  blades. 
Finally,  the  two-dimensional  grid  is  completed  by  generating  a 
distribution  of  points  in  the  0  direction  as  was  done  when  the  pseudo- 
streamsheets  were  generated.  When  completed  for  each  streamsheet,  the 
stacked  two-dimensional  grids  form  a  smooth  three-dimensional  algebraic 
periodic  H-grid  for  a  turbomachinery  blade  row.  An  example  is  given  in 
Figure  A.1.  It  should  be  noted  that  H-grids  have  a  geometric 
discontinuity  at  the  leading  edge  and  so  the  computation  in  that  region 
may  be  less  accurate.  A*so,  H-grids  are  not  recommended  for  blunt 
leading  edge  airfoils. 
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